﻿#include "HybridManufacturing.h"
#include <stdlib.h>
#include <igl/readOBJ.h>
#include <igl/writeOBJ.h>
#include <array>

bool check_inside_2(Point_2 pt, Point_2 *pgn_begin, Point_2 *pgn_end, K traits)
{
	switch (CGAL::bounded_side_2(pgn_begin, pgn_end, pt, traits))
	{
	case CGAL::ON_BOUNDED_SIDE:
		return true;
		break;
	case CGAL::ON_BOUNDARY:
		return true;
		break;
	case CGAL::ON_UNBOUNDED_SIDE:
		return false;
		break;
	}
}

HybridManufacturing::HybridManufacturing(std::string file_name, std::string suf, Eigen::MatrixXd V, Eigen::MatrixXi F, Eigen::MatrixXd NN)
{
	this->V = V;
	this->F = F;
	this->NN = NN;
	this->file_name = file_name;
	this->suf = suf;
}

HybridManufacturing::~HybridManufacturing()
{
}

void HybridManufacturing::GetVoronoiCells()
{
	int start_time = clock();
	double threshold_z = 0;
	vector<vector<Eigen::Vector3d>> all_lines;
	double min_z = MAX_D;

	for (int i = 0; i < V.rows(); i++)
		min_z = min(min_z, V(i, 2));
	for (int i = 0; i < V.rows(); i++)
	{
		if (V(i, 2) - min_z > threshold_z)
		{
			VoronoiCell new_cell;
			// new_cell.site_normal.m_x = new_cell.site_normal.m_y = new_cell.site_normal.m_z = 0;
			vector<Eigen::MatrixXd> all_boundary_V;
			vector<int> index_of_adjacent_F;
			index_of_adjacent_F.clear();
			vector<int> index_of_adjacent_V;
			index_of_adjacent_V.clear();

			//======================================独立===========================================

			// 每个顶点都可以单独地遍历所有的面，这个过程是可以并行的
			// 可以用 cuda 加速

			for (int j = 0; j < F.rows(); j++)
			{
				for (int k = 0; k < 3; k++)
				{
					if (F(j, k) == i)
					{
						index_of_adjacent_F.push_back(j);
					}
				}
			}
			//============================

			int index_begin, index_end;
			int index_first_begin;
			for (int j = 0; j < 3; j++)
				if (F(index_of_adjacent_F[0], j) == i)
				{
					index_begin = F(index_of_adjacent_F[0], (j + 1) % 3);
					index_end = F(index_of_adjacent_F[0], (j + 2) % 3);
					index_first_begin = index_begin;
					break;
				}
			index_of_adjacent_V.push_back(index_begin);
			int cont_num = 1;
			while (index_end != index_first_begin)
			{
				for (int j = cont_num; j < index_of_adjacent_F.size(); j++)
				{
					for (int k = 0; k < 3; k++)
					{
						if (F(index_of_adjacent_F[j], k) == i)
						{
							if (F(index_of_adjacent_F[j], (k + 1) % 3) == index_end)
							{
								index_begin = F(index_of_adjacent_F[j], (k + 1) % 3);
								index_end = F(index_of_adjacent_F[j], (k + 2) % 3);
								index_of_adjacent_V.push_back(index_begin);
								swap(index_of_adjacent_F[j], index_of_adjacent_F[cont_num]);
								cont_num++;
							}
							else if (F(index_of_adjacent_F[j], (k + 2) % 3) == index_end)
							{
								index_begin = F(index_of_adjacent_F[j], (k + 2) % 3);
								index_end = F(index_of_adjacent_F[j], (k + 1) % 3);
								index_of_adjacent_V.push_back(index_begin);
								swap(index_of_adjacent_F[j], index_of_adjacent_F[cont_num]);
								cont_num++;
							}
						}
					}
				}
			}
			for (int j = 0; j < index_of_adjacent_F.size(); j++)
			{
				double na = (V(F(index_of_adjacent_F[j], 1), 1) - V(F(index_of_adjacent_F[j], 0), 1)) * (V(F(index_of_adjacent_F[j], 2), 2) - V(F(index_of_adjacent_F[j], 0), 2)) - (V(F(index_of_adjacent_F[j], 1), 2) - V(F(index_of_adjacent_F[j], 0), 2)) * (V(F(index_of_adjacent_F[j], 2), 1) - V(F(index_of_adjacent_F[j], 0), 1));
				double nb = (V(F(index_of_adjacent_F[j], 1), 2) - V(F(index_of_adjacent_F[j], 0), 2)) * (V(F(index_of_adjacent_F[j], 2), 0) - V(F(index_of_adjacent_F[j], 0), 0)) - (V(F(index_of_adjacent_F[j], 1), 0) - V(F(index_of_adjacent_F[j], 0), 0)) * (V(F(index_of_adjacent_F[j], 2), 2) - V(F(index_of_adjacent_F[j], 0), 2));
				double nc = (V(F(index_of_adjacent_F[j], 1), 0) - V(F(index_of_adjacent_F[j], 0), 0)) * (V(F(index_of_adjacent_F[j], 2), 1) - V(F(index_of_adjacent_F[j], 0), 1)) - (V(F(index_of_adjacent_F[j], 1), 1) - V(F(index_of_adjacent_F[j], 0), 1)) * (V(F(index_of_adjacent_F[j], 2), 0) - V(F(index_of_adjacent_F[j], 0), 0));
				Eigen::Vector3d vn(na, nb, nc);
				// new_cell.site_normal.m_x += vn[0];
				// new_cell.site_normal.m_y += vn[1];
				// new_cell.site_normal.m_z += vn[2];
				Eigen::Matrix3d rotMatrix;

				Eigen::Vector3d vectorBefore(0, 0, 1);
				rotMatrix = Eigen::Quaterniond::FromTwoVectors(vectorBefore, vn).toRotationMatrix();
				vector<Eigen::MatrixXd> current_V;
				Eigen::MatrixXd temp_V;
				temp_V.resize(3, 1);
				temp_V(0, 0) = V(i, 0);
				temp_V(1, 0) = V(i, 1);
				temp_V(2, 0) = V(i, 2);
				current_V.push_back(temp_V);
				// faith
				if (index_of_adjacent_V.size() > j)
				{
					temp_V(0, 0) = V(index_of_adjacent_V[j], 0);
				}
				else
				{
					throw std::out_of_range("Index out of range");
				}
				temp_V(0, 0) = V(index_of_adjacent_V[j], 0);
				temp_V(1, 0) = V(index_of_adjacent_V[j], 1);
				temp_V(2, 0) = V(index_of_adjacent_V[j], 2);
				current_V.push_back(temp_V);
				temp_V(0, 0) = V(index_of_adjacent_V[(j + 1) % index_of_adjacent_V.size()], 0);
				temp_V(1, 0) = V(index_of_adjacent_V[(j + 1) % index_of_adjacent_V.size()], 1);
				temp_V(2, 0) = V(index_of_adjacent_V[(j + 1) % index_of_adjacent_V.size()], 2);
				current_V.push_back(temp_V);
				for (int k = 0; k < current_V.size(); k++)
				{
					current_V[k] = rotMatrix.inverse() * current_V[k];
				}

				/*double x1 = current_V[0](0, 0), x2 = current_V[1](0, 0);
				double y1 = current_V[0](1, 0), y2 = current_V[1](1, 0);
				double B1 = 1, A1 = (x1 - x2) / (y1 - y2); double C1 = -((y1 + y2) / 2 + A1 * (x1 + x2) / 2);
				x1 = current_V[0](0, 0), x2 = current_V[2](0, 0);
				y1 = current_V[0](1, 0), y2 = current_V[2](1, 0);
				double B2 = 1, A2 = (x1 - x2) / (y1 - y2); double C2 = -((y1 + y2) / 2 + A2 * (x1 + x2) / 2);

				Eigen::MatrixXd new_V;
				new_V.resize(3, 1);
				new_V(0, 0) = (C2 * B1 - C1 * B2) / (A1 * B2 - A2 * B1);
				new_V(1, 0) = (C1 * A2 - C2 * A1) / (A1 * B2 - A2 * B1);
				new_V(2, 0) = current_V[0](2, 0);*/

				Eigen::MatrixXd new_V;
				new_V.resize(3, 1);
				new_V(0, 0) = (current_V[0](0, 0) + current_V[1](0, 0) + current_V[2](0, 0)) / 3;
				new_V(1, 0) = (current_V[0](1, 0) + current_V[1](1, 0) + current_V[2](1, 0)) / 3;
				new_V(2, 0) = current_V[0](2, 0);

				rotMatrix = Eigen::Quaterniond::FromTwoVectors(vn, vectorBefore).toRotationMatrix();
				new_V = rotMatrix.inverse() * new_V;
				all_boundary_V.push_back(new_V);

				if (all_boundary_V.size() == 2)
				{
					Eigen::Vector3d v1 = Eigen::Vector3d(V(i, 0), V(i, 1), V(i, 2));
					Eigen::Vector3d v2 = Eigen::Vector3d(all_boundary_V[0](0, 0), all_boundary_V[0](1, 0), all_boundary_V[0](2, 0));
					Eigen::Vector3d v3 = Eigen::Vector3d(all_boundary_V[1](0, 0), all_boundary_V[1](1, 0), all_boundary_V[1](2, 0));
					// double ans = (v2.x() - v1.x()) * (v2.y() - v3.y()) - (v2.y() - v1.y()) * (v2.x() - v3.x());
					// if (ans > 0)	//is clockwise
					//	swap(v2, v3);
					double na = (v2.y() - v1.y()) * (v3.z() - v1.z()) - (v2.z() - v1.z()) * (v3.y() - v1.y());
					double nb = (v2.z() - v1.z()) * (v3.x() - v1.x()) - (v2.x() - v1.x()) * (v3.z() - v1.z());
					double nc = (v2.x() - v1.x()) * (v3.y() - v1.y()) - (v2.y() - v1.y()) * (v3.x() - v1.x());
					Eigen::Vector3d vn_2(na, nb, nc);
					if (vn.dot(vn_2) < 0)
						swap(all_boundary_V[0], all_boundary_V[1]);
				}
			}

			// new_cell.site_normal /= index_of_adjacent_F.size();
			// new_cell.site_normal.Normalized();
			// new_cell.index_of_V = i;
			new_cell.is_available = true;
			new_cell.site = i;
			new_cell.adjacent_cells = index_of_adjacent_V;
			vector<Eigen::Vector3d> temp_lines;
			all_lines.push_back(temp_lines);
			for (int j = 0; j < all_boundary_V.size(); j++)
			{
				Eigen::Vector3d temp_vec(all_boundary_V[j](0, 0), all_boundary_V[j](1, 0), all_boundary_V[j](2, 0));
				all_lines[all_lines.size() - 1].push_back(temp_vec);
				new_cell.all_points_in_polygon.push_back(Vector3(all_boundary_V[j](0, 0), all_boundary_V[j](1, 0), all_boundary_V[j](2, 0)));
			}
			all_voronoi_cells.push_back(new_cell);
		}
		else
		{
			V_bottom.push_back(Vector3(V(i, 0), V(i, 1), V(i, 2)));
			VoronoiCell new_cell;
			new_cell.is_available = false;
			all_voronoi_cells.push_back(new_cell);
		}
	}
	int end_time = clock();

	Visual vis;
	if (open_vis_voronoi == true)
		vis.generateModelForRendering_8(all_lines, file_name + "_voronoi.obj");
}

bool HybridManufacturing::JudgeIntersectTriangleAndCircle(Vector3 v1, Vector3 v2, Vector3 v3, Vector3 center, double r)
{
	return false;
}

int HybridManufacturing::CollisionDetectionForSubtractiveManufacturing(cutter the_nozzle, clock_t start_time_collision)
{
	ofstream cutter_file(file_name + "_cutter.txt");

	if (if_keep_other_info)
		cutter_file << the_nozzle.cylinder_r << " " << the_nozzle.cylinder_height << " " << the_nozzle.carriage_r << " " << the_nozzle.carriage_height << endl;
	clock_t start_time, end_time;
	clock_t start_time_2, end_time_2;
	clock_t start_time_3, end_time_3;
	clock_t start_time_4, end_time_4;
	double sum_time = 0;
	int test_cont_num = 0;

	start_time = clock();

	string file_name_2 = file_name;
	int cont_num = 0;
	sampling_subtractive.OrientationSamplePoints();
	Eigen::Matrix3d rotMatrix;
	std::vector<bool> flag_accessible_points;
	std::vector<bool> flag_covering_points;

	ofstream file_out(file_name + "_collision_detection.txt");

	vector<bool> is_accessible_points(V.rows());
	for (int i = 0; i < V.rows(); i++)
		is_accessible_points[i] = true;

	flag_accessible_points.resize(all_voronoi_cells.size());
	flag_covering_points.resize(all_voronoi_cells.size());
	flag_voronoi_has_been_printed.resize(all_voronoi_cells.size());

	for (int i = 0; i < all_voronoi_cells.size(); i++)
	{
		flag_accessible_points[i] = false;
		flag_voronoi_has_been_printed[i] = true;
		flag_covering_points[i] = false;
	}

	for (int ori = 0; ori < sampling_subtractive.sample_points.size(); ori++)
	{

		///////////////////rotate/////////////////////
		std::vector<Eigen::MatrixXd> temp_V;
		temp_V.resize(V.rows());
		for (int i = 0; i < V.rows(); i++)
		{
			temp_V[i].resize(3, 1);
			temp_V[i](0, 0) = V.row(i).x();
			temp_V[i](1, 0) = V.row(i).y();
			temp_V[i](2, 0) = V.row(i).z();
		}
		vector<std::vector<Eigen::MatrixXd>> temp_new_V;
		temp_new_V.resize(V.rows());
		for (int i = 0; i < V.rows(); i++)
		{
			temp_new_V[i].resize(all_voronoi_cells[i].all_points_in_polygon.size());
			for (int k = 0; k < temp_new_V[i].size(); k++)
			{
				temp_new_V[i][k].resize(3, 1);
				temp_new_V[i][k](0, 0) = all_voronoi_cells[i].all_points_in_polygon[k].m_x;
				temp_new_V[i][k](1, 0) = all_voronoi_cells[i].all_points_in_polygon[k].m_y;
				temp_new_V[i][k](2, 0) = all_voronoi_cells[i].all_points_in_polygon[k].m_z;
			}
		}

		Eigen::Vector3d vectorBefore(0, 0, 1);
		Eigen::Vector3d vectorAfter(sampling_subtractive.sample_points[ori].x, sampling_subtractive.sample_points[ori].y, sampling_subtractive.sample_points[ori].z);
		rotMatrix = Eigen::Quaterniond::FromTwoVectors(vectorBefore, vectorAfter).toRotationMatrix();
		for (int i = 0; i < V.rows(); i++)
			temp_V[i] = rotMatrix.inverse() * temp_V[i];
		for (int i = 0; i < V.rows(); i++)
			for (int j = 0; j < temp_new_V[i].size(); j++)
				temp_new_V[i][j] = rotMatrix.inverse() * temp_new_V[i][j];
		//////////////////////////////////////////////

		/////////////////calculate normal//////////////////
		// vector<vector<Vector3>> all_normal_of_triangles_in_cells;
		vector<Vector3> all_normal_of_cells;
		// all_normal_of_triangles_in_cells.resize(all_voronoi_cells.size());
		all_normal_of_cells.resize(all_voronoi_cells.size());

		vector<vector<cv::Point3d>> temp_vis;
		for (int i = 0; i < all_voronoi_cells.size(); i++)
		{
			if (all_voronoi_cells[i].is_available == true)
			{
				// all_normal_of_triangles_in_cells[i].resize(all_voronoi_cells[i].all_points_in_polygon.size());
				all_normal_of_cells[i].m_x = all_normal_of_cells[i].m_y = all_normal_of_cells[i].m_z = 0;
				for (int j = 0; j < 1; j++)
				{
					Vector3 v1 = Vector3(temp_V[all_voronoi_cells[i].site](0, 0), temp_V[all_voronoi_cells[i].site](1, 0), temp_V[all_voronoi_cells[i].site](2, 0));
					Vector3 v2 = Vector3(temp_new_V[i][j](0, 0), temp_new_V[i][j](1, 0), temp_new_V[i][j](2, 0));
					Vector3 v3 = Vector3(temp_new_V[i][(j + 1) % all_voronoi_cells[i].all_points_in_polygon.size()](0, 0), temp_new_V[i][(j + 1) % all_voronoi_cells[i].all_points_in_polygon.size()](1, 0), temp_new_V[i][(j + 1) % all_voronoi_cells[i].all_points_in_polygon.size()](2, 0));
					double ans = (v2.m_x - v1.m_x) * (v2.m_y - v3.m_y) - (v2.m_y - v1.m_y) * (v2.m_x - v3.m_x);
					// if (ans > 0)	//is clockwise
					//	swap(v2, v3);
					double na = (v2.m_y - v1.m_y) * (v3.m_z - v1.m_z) - (v2.m_z - v1.m_z) * (v3.m_y - v1.m_y);
					double nb = (v2.m_z - v1.m_z) * (v3.m_x - v1.m_x) - (v2.m_x - v1.m_x) * (v3.m_z - v1.m_z);
					double nc = (v2.m_x - v1.m_x) * (v3.m_y - v1.m_y) - (v2.m_y - v1.m_y) * (v3.m_x - v1.m_x);
					Vector3 vn(na, nb, nc);
					// all_normal_of_triangles_in_cells[i][j] = vn;
					all_normal_of_cells[i].m_x += vn.m_x;
					all_normal_of_cells[i].m_y += vn.m_y;
					all_normal_of_cells[i].m_z += vn.m_z;
				}
				all_normal_of_cells[i] /= all_voronoi_cells[i].all_points_in_polygon.size();
				all_normal_of_cells[i].Normalized();

				all_voronoi_cells[i].all_normal_in_all_ori.push_back(all_normal_of_cells[i]);

				vector<cv::Point3d> temp_vec;
				temp_vis.push_back(temp_vec);
				cv::Point3d v_site = cv::Point3d(temp_V[all_voronoi_cells[i].site](0, 0), temp_V[all_voronoi_cells[i].site](1, 0), temp_V[all_voronoi_cells[i].site](2, 0));
				cv::Point3d v_normal = cv::Point3d(all_normal_of_cells[i].m_x * 3 + v_site.x, all_normal_of_cells[i].m_y * 3 + v_site.y, all_normal_of_cells[i].m_z * 3 + v_site.z);
				temp_vis[temp_vis.size() - 1].push_back(v_site);
				temp_vis[temp_vis.size() - 1].push_back(v_normal);
			}
		}
		if (ori == 0)
		{
			Visual vis_normal;
			// vis_normal.generateModelForRendering_9(temp_vis, file_name);
		}
		//////////////////////////////////////////////////

		///////////////////collision detection////////////////////////
		int cont_accessible_points = 0, cont_unaccessible_points = 0;
		double circle_r;

		vector<double> max_z_of_cells(all_voronoi_cells.size());

		for (int i = 0; i < all_voronoi_cells.size(); i++)
		{
			max_z_of_cells[i] = MIN_D;
			for (int j = 0; j < all_voronoi_cells[i].all_points_in_polygon.size(); j++)
				max_z_of_cells[i] = max(max_z_of_cells[i], temp_new_V[i][j](2, 0));
		}

		for (int i = 0; i < all_voronoi_cells.size(); i++)
		{
			if (flag_accessible_points[i] == true)
			{
				cont_accessible_points++;
				continue;
			}
			cont_unaccessible_points++;
			bool jud_collision = false;
			Vector3 center_point;
			if (all_voronoi_cells[i].is_available == true)
			{
				center_point.m_x = temp_V[all_voronoi_cells[i].site](0, 0) + the_nozzle.cylinder_r * all_normal_of_cells[i].m_x;
				center_point.m_y = temp_V[all_voronoi_cells[i].site](1, 0) + the_nozzle.cylinder_r * all_normal_of_cells[i].m_y;
				center_point.m_z = temp_V[all_voronoi_cells[i].site](2, 0) + the_nozzle.cylinder_r * all_normal_of_cells[i].m_z;
			}

			for (int ii = 0; ii < all_voronoi_cells.size(); ii++)
			{
				if (i != ii && all_voronoi_cells[i].is_available == true && all_voronoi_cells[ii].is_available == true)
				{

					// 碰撞检测，三角形与竖直圆柱，暂时用顶点判断//
					bool jud_collision_2 = false;
					// Vector3 v1 = Vector3 (temp_V[all_voronoi_cells[ii].site](0,0), temp_V[all_voronoi_cells[ii].site](1, 0), temp_V[all_voronoi_cells[ii].site](2, 0));
					Vector3 v_boundary;
					if (max_z_of_cells[ii] - center_point.m_z <= 0)
					{
						jud_collision = false;
						continue;
					}
					else if (max_z_of_cells[ii] - center_point.m_z > the_nozzle.cylinder_height + the_nozzle.ball_r + the_nozzle.carriage_height)
					{
						jud_collision = true;
						break;
					}
					else if (max_z_of_cells[ii] - center_point.m_z > the_nozzle.cylinder_height + the_nozzle.ball_r)
					{
						if (sqrt(pow(temp_new_V[ii][0](0, 0) - center_point.m_x, 2) + pow(temp_new_V[ii][0](1, 0) - center_point.m_y, 2)) > the_nozzle.carriage_r + 5)
						{
							jud_collision = false;
							continue;
						}
					}
					else if (max_z_of_cells[ii] - center_point.m_z <= the_nozzle.cylinder_height + the_nozzle.ball_r)
					{
						if (sqrt(pow(temp_new_V[ii][0](0, 0) - center_point.m_x, 2) + pow(temp_new_V[ii][0](1, 0) - center_point.m_y, 2)) > the_nozzle.cylinder_r + 5)
						{
							jud_collision = false;
							continue;
						}
					}
					for (int j = 0; j < all_voronoi_cells[ii].all_points_in_polygon.size(); j++)
					{
						v_boundary = Vector3(temp_new_V[ii][j](0, 0), temp_new_V[ii][j](1, 0), temp_new_V[ii][j](2, 0));
						double diff_z = v_boundary.m_z - center_point.m_z;
						if (diff_z <= 0)
						{
							continue;
						}
						else if (diff_z <= the_nozzle.cylinder_r + the_nozzle.cylinder_height)
						{
							if (pow(v_boundary.m_x - center_point.m_x, 2) + pow(v_boundary.m_y - center_point.m_y, 2) - pow(the_nozzle.cylinder_r, 2) < 0)
							{
								jud_collision_2 = true;
								break;
							}
						}
						else if (diff_z <= the_nozzle.cylinder_r + the_nozzle.cylinder_height + the_nozzle.carriage_height)
						{
							if (pow(v_boundary.m_x - center_point.m_x, 2) + pow(v_boundary.m_y - center_point.m_y, 2) - pow(the_nozzle.carriage_r, 2) < 0)
							{
								jud_collision_2 = true;
								break;
							}
						}
						// v3 = Vector3(temp_new_V[ii][(j + 1) % all_voronoi_cells[ii].all_points_in_polygon.size()](0, 0), temp_new_V[ii][(j + 1) % all_voronoi_cells[ii].all_points_in_polygon.size()](1, 0), temp_new_V[ii][(j + 1) % all_voronoi_cells[ii].all_points_in_polygon.size()](2, 0));
						// jud_collision_2 = JudgeIntersectTriangleAndCircle(v1, v2 , v3, center_point, the_nozzle.cylinder_r);
					}
					if (jud_collision_2 == true)
					{
						jud_collision = true;
						break;
					}
					//////////////////
				}
			}
			if (jud_collision == false)
			{
				flag_accessible_points[i] = true;
			}
		}
		num_inaccessible_points = cont_unaccessible_points;
		if (if_debug == true)
		{
			std::cout << "id of orientation:" << ori << endl;
			std::cout << "number of accessible points:" << cont_accessible_points << endl;
			std::cout << "number of unaccessible points:" << cont_unaccessible_points << endl
					  << endl;
		}
		// break;
	}
	end_time = clock();

	//*******************************************//
	/////////////////find area S//////////////////
	//*******************************************//
	/*start_time_2 = clock();
	int cont_number_2 = 0;
	std::vector<Eigen::MatrixXd> vis_red_points;
	std::vector<Eigen::MatrixXd> vis_green_points;
	std::vector<Eigen::MatrixXd> temp_V_vis;
	std::vector<Eigen::MatrixXd> temp_V;
	temp_V.resize(V.rows());
	temp_V_vis.resize(V.rows());
	for (int i = 0; i < V.rows(); i++) {
		temp_V[i].resize(3, 1);
		temp_V[i](0, 0) = V.row(i).x();
		temp_V[i](1, 0) = V.row(i).y();
		temp_V[i](2, 0) = V.row(i).z();
	}
	temp_V_vis = temp_V;
	for (int i = 0; i < all_voronoi_cells.size(); i++) {
		if (flag_accessible_points[i] == true) {
			continue;
		}
		vis_red_points.push_back(temp_V_vis[i]);
		is_accessible_points[i] = false;
	}
	if (open_vis_red_points == true)
		creat_red_ball(file_name, vis_red_points);
	for (int i = 0; i < V.rows(); i++) {
		if (is_accessible_points[i] == true) {
			file_out << V.row(i).x() << " " << V.row(i).y() << " " << V.row(i).z() << " " << NN.row(i).x() << " " << NN.row(i).y() << " " << NN.row(i).z() << " 1" << endl;
		}
		else {
			file_out << V.row(i).x() << " " << V.row(i).y() << " " << V.row(i).z() << " " << NN.row(i).x() << " " << NN.row(i).y() << " " << NN.row(i).z() << " 0" << endl;
		}
	}*/

	clock_t end_time_collision = clock();

	clock_t start_time_obscure = clock();

	// begin

	int cont_number_2 = 0;
	std::vector<Eigen::MatrixXd> vis_red_points;
	std::vector<Eigen::MatrixXd> vis_green_points;
	std::vector<Eigen::MatrixXd> temp_V_vis;
	std::vector<Eigen::MatrixXd> vis_most_green_points;
	std::vector<int> new_flag_covering_points(V.rows(), 0);
	vector<double> color_map;
	temp_V_vis.resize(V.rows());

	if (open_calculate_collision_points == true)
	{
		sampling_subtractive.OrientationSamplePoints_2();
		vector<int> cont_covered_points(V.rows());

		for (int i = 0; i < V.rows(); i++)
			cont_covered_points[i] = 0;
		int cont_all_covered_points = 0;
		for (int ori = 0; ori < sampling_subtractive.sample_points.size(); ori++)
		{
			/*cout << ori << endl;*/
			std::vector<Eigen::MatrixXd> temp_V;
			temp_V.resize(V.rows());
			for (int i = 0; i < V.rows(); i++)
			{
				temp_V[i].resize(3, 1);
				temp_V[i](0, 0) = V.row(i).x();
				temp_V[i](1, 0) = V.row(i).y();
				temp_V[i](2, 0) = V.row(i).z();
			}
			vector<std::vector<Eigen::MatrixXd>> temp_new_V;
			temp_new_V.resize(V.rows());
			for (int i = 0; i < V.rows(); i++)
			{
				temp_new_V[i].resize(all_voronoi_cells[i].all_points_in_polygon.size());
				for (int k = 0; k < temp_new_V[i].size(); k++)
				{
					temp_new_V[i][k].resize(3, 1);
					temp_new_V[i][k](0, 0) = all_voronoi_cells[i].all_points_in_polygon[k].m_x;
					temp_new_V[i][k](1, 0) = all_voronoi_cells[i].all_points_in_polygon[k].m_y;
					temp_new_V[i][k](2, 0) = all_voronoi_cells[i].all_points_in_polygon[k].m_z;
				}
			}
			temp_V_vis = temp_V;

			Eigen::Vector3d vectorBefore(0, 0, 1);
			Eigen::Vector3d vectorAfter(sampling_subtractive.sample_points[ori].x, sampling_subtractive.sample_points[ori].y, sampling_subtractive.sample_points[ori].z);
			rotMatrix = Eigen::Quaterniond::FromTwoVectors(vectorBefore, vectorAfter).toRotationMatrix();
			for (int i = 0; i < V.rows(); i++)
				temp_V[i] = rotMatrix.inverse() * temp_V[i];
			for (int i = 0; i < V.rows(); i++)
				for (int j = 0; j < temp_new_V[i].size(); j++)
					temp_new_V[i][j] = rotMatrix.inverse() * temp_new_V[i][j];

			// cout << "b" << endl;
			/////////////////calculate normal//////////////////
			// vector<vector<Vector3>> all_normal_of_triangles_in_cells;
			vector<Vector3> all_normal_of_cells;
			// all_normal_of_triangles_in_cells.resize(all_voronoi_cells.size());
			all_normal_of_cells.resize(all_voronoi_cells.size());
			for (int i = 0; i < all_voronoi_cells.size(); i++)
			{
				if (all_voronoi_cells[i].is_available == true && flag_accessible_points[i] == false)
				{
					// all_normal_of_triangles_in_cells[i].resize(all_voronoi_cells[i].all_points_in_polygon.size());
					all_normal_of_cells[i].m_x = all_normal_of_cells[i].m_y = all_normal_of_cells[i].m_z = 0;
					for (int j = 0; j < 1; j++)
					{
						Vector3 v1 = Vector3(temp_V[all_voronoi_cells[i].site](0, 0), temp_V[all_voronoi_cells[i].site](1, 0), temp_V[all_voronoi_cells[i].site](2, 0));
						Vector3 v2 = Vector3(temp_new_V[i][j](0, 0), temp_new_V[i][j](1, 0), temp_new_V[i][j](2, 0));
						Vector3 v3 = Vector3(temp_new_V[i][(j + 1) % all_voronoi_cells[i].all_points_in_polygon.size()](0, 0), temp_new_V[i][(j + 1) % all_voronoi_cells[i].all_points_in_polygon.size()](1, 0), temp_new_V[i][(j + 1) % all_voronoi_cells[i].all_points_in_polygon.size()](2, 0));
						double ans = (v2.m_x - v1.m_x) * (v2.m_y - v3.m_y) - (v2.m_y - v1.m_y) * (v2.m_x - v3.m_x);
						// if (ans > 0)	//is clockwise
						//	swap(v2, v3);
						double na = (v2.m_y - v1.m_y) * (v3.m_z - v1.m_z) - (v2.m_z - v1.m_z) * (v3.m_y - v1.m_y);
						double nb = (v2.m_z - v1.m_z) * (v3.m_x - v1.m_x) - (v2.m_x - v1.m_x) * (v3.m_z - v1.m_z);
						double nc = (v2.m_x - v1.m_x) * (v3.m_y - v1.m_y) - (v2.m_y - v1.m_y) * (v3.m_x - v1.m_x);
						Vector3 vn(na, nb, nc);
						// all_normal_of_triangles_in_cells[i][j] = vn;
						all_normal_of_cells[i].m_x += vn.m_x;
						all_normal_of_cells[i].m_y += vn.m_y;
						all_normal_of_cells[i].m_z += vn.m_z;
					}
					all_normal_of_cells[i] /= all_voronoi_cells[i].all_points_in_polygon.size();
					all_normal_of_cells[i].Normalized();
				}
			}
			//////////////////////////////////////////////////
			// cout << "c" << endl;
			// ofstream vis_unaccessible_points(file_name + "_unaccessible_points.obj");

			file_name_2 = file_name + to_string(cont_num);
			vector<double> max_z_of_cells(all_voronoi_cells.size());
			for (int i = 0; i < all_voronoi_cells.size(); i++)
			{
				max_z_of_cells[i] = MIN_D;
				for (int j = 0; j < all_voronoi_cells[i].all_points_in_polygon.size(); j++)
					max_z_of_cells[i] = max(max_z_of_cells[i], temp_new_V[i][j](2, 0));
			}

			double circle_r;
			int cont_unaccessible_points = 0;
			for (int i = 0; i < all_voronoi_cells.size(); i++)
			{
				if (flag_accessible_points[i] == true)
				{
					continue;
				}

				Vector3 center_point;
				if (all_voronoi_cells[i].is_available == true)
				{
					center_point.m_x = temp_V[all_voronoi_cells[i].site](0, 0) + the_nozzle.cylinder_r * all_normal_of_cells[i].m_x;
					center_point.m_y = temp_V[all_voronoi_cells[i].site](1, 0) + the_nozzle.cylinder_r * all_normal_of_cells[i].m_y;
					center_point.m_z = temp_V[all_voronoi_cells[i].site](2, 0) + the_nozzle.cylinder_r * all_normal_of_cells[i].m_z;
				}

				// vis_unaccessible_points << "v " << temp_V[i](0,0) << " " << temp_V[i](1, 0) << " " << temp_V[i](2, 0) << endl;
				start_time_4 = clock();
				if (ori == 0)
				{
					vis_red_points.push_back(temp_V_vis[i]);
					is_accessible_points[i] = false;
				}
				end_time_4 = clock();
				sum_time += (double(end_time_4) - double(start_time_4)) / CLOCKS_PER_SEC;
				for (int ii = 0; ii < all_voronoi_cells.size(); ii++)
				{
					if (i != ii && all_voronoi_cells[i].is_available == true && all_voronoi_cells[ii].is_available == true)
					{
						bool jud_collision_2 = false;
						// Vector3 v1 = Vector3 (temp_V[all_voronoi_cells[ii].site](0,0), temp_V[all_voronoi_cells[ii].site](1, 0), temp_V[all_voronoi_cells[ii].site](2, 0));
						Vector3 v_boundary;
						if (max_z_of_cells[ii] - center_point.m_z <= 0)
						{
							continue;
						}
						if (max_z_of_cells[ii] - center_point.m_z > the_nozzle.cylinder_height + the_nozzle.ball_r + the_nozzle.carriage_height)
						{
							area_S temp_area_S(ii, ori);
							cont_covered_points[ii]++;
							cont_all_covered_points++;
							continue;
						}
						else if (max_z_of_cells[ii] - center_point.m_z > the_nozzle.cylinder_height + the_nozzle.ball_r)
						{
							if (sqrt(pow(temp_new_V[ii][0](0, 0) - center_point.m_x, 2) + pow(temp_new_V[ii][0](1, 0) - center_point.m_y, 2)) > the_nozzle.carriage_r + 5)
							{
								continue;
							}
						}
						else if (max_z_of_cells[ii] - center_point.m_z <= the_nozzle.cylinder_height + the_nozzle.ball_r)
						{
							if (sqrt(pow(temp_new_V[ii][0](0, 0) - center_point.m_x, 2) + pow(temp_new_V[ii][0](1, 0) - center_point.m_y, 2)) > the_nozzle.cylinder_r + 5)
							{
								continue;
							}
						}
						for (int j = 0; j < all_voronoi_cells[ii].all_points_in_polygon.size(); j++)
						{
							v_boundary = Vector3(temp_new_V[ii][j](0, 0), temp_new_V[ii][j](1, 0), temp_new_V[ii][j](2, 0));
							double diff_z = v_boundary.m_z - center_point.m_z;
							if (diff_z <= 0)
								continue;
							else if (diff_z <= the_nozzle.cylinder_r + the_nozzle.cylinder_height)
							{
								if (pow(v_boundary.m_x - center_point.m_x, 2) + pow(v_boundary.m_y - center_point.m_y, 2) - pow(the_nozzle.cylinder_r, 2) < 0)
								{
									jud_collision_2 = true;
									break;
								}
							}
							else if (diff_z <= the_nozzle.cylinder_r + the_nozzle.cylinder_height + the_nozzle.carriage_height)
							{
								if (pow(v_boundary.m_x - center_point.m_x, 2) + pow(v_boundary.m_y - center_point.m_y, 2) - pow(the_nozzle.carriage_r, 2) < 0)
								{
									jud_collision_2 = true;
									break;
								}
							}
							// v3 = Vector3(temp_new_V[ii][(j + 1) % all_voronoi_cells[ii].all_points_in_polygon.size()](0, 0), temp_new_V[ii][(j + 1) % all_voronoi_cells[ii].all_points_in_polygon.size()](1, 0), temp_new_V[ii][(j + 1) % all_voronoi_cells[ii].all_points_in_polygon.size()](2, 0));
							// jud_collision_2 = JudgeIntersectTriangleAndCircle(v1, v2 , v3, center_point, the_nozzle.cylinder_r);
						}
						if (jud_collision_2 == true)
						{
							cont_covered_points[ii]++;
							cont_all_covered_points++;
						}
						/////////////////////////////////////////////
					}
				}
			}

			// cout << "d" << endl;
			Visual Vis;
			// Vis.generateModelForRendering(temp_layers, file_name_2);
			// Vis.generateModelForRendering_3(vectorAfter, file_name_2, vis_points);
			// break;

			cont_num++;
		}
		//*******************************************//
		// cout << endl << all_the_area_S.size() << endl;
		// cout << endl << cont_number_2 << endl;
		end_time_2 = clock();

		////////////get covering points(green points)////////////
		// only consider covering for area S,save the index of all_the_area_S
		start_time_3 = clock();
		cout << test_cont_num << "&&&&" << endl;

		std::vector<std::pair<int, int>> size_with_index(V.rows());
		// 计算每个点集的大小和索引
		for (int i = 0; i < V.rows(); i++)
		{
			size_with_index[i] = {cont_covered_points[i], i};
		}
		// 按照点集大小从大到小排序
		std::sort(size_with_index.begin(), size_with_index.end(), std::greater<>());
		// 计算前10%的数量
		int top_10_percent_count = 0;
		for (int i = 0; i < V.rows(); i++)
		{
			if (cont_covered_points[i] != 0)
				top_10_percent_count++;
		}
		top_10_percent_count = int(top_10_percent_count / 10);
		// 将前10%对应的 flag_covering_points 设为1
		if (cont_all_covered_points != 0)
			for (int k = 0; k < top_10_percent_count; k++)
			{
				int index = size_with_index[k].second;				// 获取原始索引
				new_flag_covering_points[index] = 1;				// 设置对应的 flag
				vis_most_green_points.push_back(temp_V_vis[index]); // 添加到可视化点集
			}

		color_map.resize(vis_most_green_points.size());
		for (int i = 0; i < vis_most_green_points.size(); i++)
			color_map[i] = 1;

		end_time_3 = clock();
	}
	else
	{
		for (int i = 0; i < V.rows(); i++)
		{
			temp_V_vis[i].resize(3, 1);
			temp_V_vis[i](0, 0) = V.row(i).x();
			temp_V_vis[i](1, 0) = V.row(i).y();
			temp_V_vis[i](2, 0) = V.row(i).z();
		}
		for (int i = 0; i < all_voronoi_cells.size(); i++)
		{
			if (flag_accessible_points[i] != true)
			{
				vis_red_points.push_back(temp_V_vis[i]);
				is_accessible_points[i] = false;
			}
		}
	}

	clock_t end_time_obscure = clock();

	/////////////////////////////////////////////////////////

	if (open_vis_red_points == true)
		creat_red_ball(file_name, vis_red_points);
	if (open_vis_green_points == true)
		creat_green_ball(file_name, vis_most_green_points, color_map);

	/////////////////////////////////////////////////////////////////////////////////////////////////
	if (open_count_red_points == true)
	{
		ofstream file_unaccessible_out(file_name + "_ori_" +suf + "_unaccessible_point_count.txt");
		int total_inaccessible_points = 0;
		int total_point_number = flag_accessible_points.size();

		for (int i = 0; i < total_point_number; i++)
		{
			if (!flag_accessible_points[i])
			{
				total_inaccessible_points++;
			}
		}
		file_unaccessible_out << total_point_number << " " << total_inaccessible_points << endl;
		double ratio = 0.0;
		if (total_point_number > 0)
		{
			ratio = double(total_inaccessible_points) / total_point_number * 100.0;
		}
		file_unaccessible_out << std::fixed << std::setprecision(4) << ratio << "%" << endl;
	}
	if (open_calculate_collision_points == true)
	{
		if (open_keep_points)
		{
			for (int i = 0; i < V.rows(); i++)
			{
				if (is_accessible_points[i] == true)
				{
					file_out << V.row(i).x() << " " << V.row(i).y() << " " << V.row(i).z() << " " << NN.row(i).x() << " " << NN.row(i).y() << " " << NN.row(i).z() << " 1 " << 0 << endl;
				}
				else
				{
					file_out << V.row(i).x() << " " << V.row(i).y() << " " << V.row(i).z() << " " << NN.row(i).x() << " " << NN.row(i).y() << " " << NN.row(i).z() << " 0 " << 0 << endl;
				}
			}
		}
	}
	else
	{
		for (int i = 0; i < V.rows(); i++)
		{
			if (is_accessible_points[i] == true)
			{
				file_out << V.row(i).x() << " " << V.row(i).y() << " " << V.row(i).z() << " " << NN.row(i).x() << " " << NN.row(i).y() << " " << NN.row(i).z() << " 1 " << endl;
			}
			else
			{
				file_out << V.row(i).x() << " " << V.row(i).y() << " " << V.row(i).z() << " " << NN.row(i).x() << " " << NN.row(i).y() << " " << NN.row(i).z() << " 0 " << endl;
			}
		}
	}

	/*std::cout << "&&&time&&& Collision detection: " << double(end_time - start_time) / CLOCKS_PER_SEC << std::endl;
	std::cout << "&&&time&&& Subtractive dependency graph: " << double(end_time_2 - start_time_2) / CLOCKS_PER_SEC << std::endl;
	std::cout << "&&&time&&& Green points generation: " << double(end_time_3 - start_time_3) / CLOCKS_PER_SEC << std::endl;*/
	/*time_build_subtractive_graph = double(end_time - start_time) / CLOCKS_PER_SEC + double(end_time_2 - start_time_2) / CLOCKS_PER_SEC + double(end_time_3 - start_time_3) / CLOCKS_PER_SEC;
	std::cout << "&&&time&&& Build subtractive graph: " << double(end_time - start_time) / CLOCKS_PER_SEC + double(end_time_2 - start_time_2) / CLOCKS_PER_SEC + double(end_time_3 - start_time_3) / CLOCKS_PER_SEC << std::endl;*/

	/*if (all_the_area_S.size() == 0) {
		cout << file_name << " 无不可达点！！！" << endl;

		Eigen::MatrixXd V;
		Eigen::MatrixXi F;
		Eigen::MatrixXi N;

		igl::readOBJ(file_name + ".obj", V, F);
		igl::writeOBJ(file_name +  "-1_0_current.obj", V, F);
		return -1;
	}*/
	return 0;
}

void HybridManufacturing::creat_red_ball(string file_name, std::vector<Eigen::MatrixXd> vis_points)
{
	Eigen::MatrixXd V_2;
	Eigen::MatrixXi F_2;

	igl::readOBJ("C:\\Users\\honor\\Desktop\\orientation\\ball_low.obj", V_2, F_2);
	ofstream all_balls(file_name + "_unaccessivle_points.obj");
	for (int i = 0; i < vis_points.size(); i++)
	{
		for (int j = 0; j < V_2.rows(); j++)
			all_balls << "v " << V_2(j, 0) + vis_points[i](0, 0) << " " << V_2(j, 1) + vis_points[i](1, 0) << " " << V_2(j, 2) + vis_points[i](2, 0) << " 0.9" << " 0.05" << " 0.05" << endl;
		for (int j = 0; j < F_2.rows(); j++)
			all_balls << "f " << F_2(j, 0) + i * V_2.rows() + 1 << " " << F_2(j, 1) + i * V_2.rows() + 1 << " " << F_2(j, 2) + i * V_2.rows() + 1 << endl;
	}
	all_balls.close();
}

void HybridManufacturing::creat_green_ball(string file_name, std::vector<Eigen::MatrixXd> vis_points, vector<double> color_map)
{
	Eigen::MatrixXd V_2;
	Eigen::MatrixXi F_2;

	igl::readOBJ("ball.obj", V_2, F_2);
	ofstream all_balls(file_name + "_covering_points.obj");
	for (int i = 0; i < vis_points.size(); i++)
	{
		for (int j = 0; j < V_2.rows(); j++)
			all_balls << "v " << V_2(j, 0) + vis_points[i](0, 0) << " " << V_2(j, 1) + vis_points[i](1, 0) << " " << V_2(j, 2) + vis_points[i](2, 0) << " " << "0.05" << " " << color_map[i] * 1.0 << " " << "0.05" << endl;
		for (int j = 0; j < F_2.rows(); j++)
			all_balls << "f " << F_2(j, 0) + i * V_2.rows() + 1 << " " << F_2(j, 1) + i * V_2.rows() + 1 << " " << F_2(j, 2) + i * V_2.rows() + 1 << endl;
	}
	all_balls.close();
}

void HybridManufacturing::SortEdges(vector<vector<area_S>> ori_all_the_area_S)
{
	for (int i = 0; i < ori_all_the_area_S.size(); i++)
	{
		for (int j = 0; j < ori_all_the_area_S[i].size(); j++)
		{
			for (int k = j + 1; k < ori_all_the_area_S[i].size(); k++)
			{
				if (ori_all_the_area_S[i][j].id_ori > ori_all_the_area_S[i][k].id_ori)
					swap(ori_all_the_area_S[i][j], ori_all_the_area_S[i][k]);
			}
		}
	}
}

double HybridManufacturing::distance_3d(cv::Point3d a, cv::Point3d b)
{
	double distance = sqrt(pow(a.x - b.x, 2) + pow(a.y - b.y, 2) + pow(a.z - b.z, 2));
	return distance;
}

double HybridManufacturing::distance_2d(cv::Point2d a, cv::Point2d b)
{
	double distance = sqrt(pow(a.x - b.x, 2) + pow(a.y - b.y, 2));
	return distance;
}

double HybridManufacturing::distance_3d_vector(VEctor a, VEctor b)
{
	double distance = sqrt(pow(a[0] - b[0], 2) + pow(a[1] - b[1], 2) + pow(a[2] - b[2], 2));
	return distance;
}

Vector3 HybridManufacturing::calculate_normal(VEctor v1, VEctor v2, VEctor v3)
{
	// double ans = (v2[0] - v1[0]) * (v2[1] - v3[1]) - (v2[1] - v1[1]) * (v2[0] - v3[0]);
	// if (ans > 0)	//is clockwise
	//	swap(v2, v3);

	double na = (v2[1] - v1[1]) * (v3[2] - v1[2]) - (v2[2] - v1[2]) * (v3[1] - v1[1]);
	double nb = (v2[2] - v1[2]) * (v3[0] - v1[0]) - (v2[0] - v1[0]) * (v3[2] - v1[2]);
	double nc = (v2[0] - v1[0]) * (v3[1] - v1[1]) - (v2[1] - v1[1]) * (v3[0] - v1[0]);

	Vector3 vn(na, nb, nc);
	return vn;
}

double HybridManufacturing::calculate_area(VEctor v1, VEctor v2, VEctor v3)
{
	Eigen::Vector3d vec1, vec2;
	vec1.x() = v2[0] - v1[0];
	vec1.y() = v2[1] - v1[1];
	vec1.z() = v2[2] - v1[2];
	vec2.x() = v3[0] - v1[0];
	vec2.y() = v3[1] - v1[1];
	vec2.z() = v3[2] - v1[2];

	double area = vec1.cross(vec2).norm() / 2;
	return area;
}

double HybridManufacturing::calculate_manufacturing_require(Slicer_2 slicer, int index, double min_z, double max_z)
{
	double average_z = (slicer.positions[slicer.triangles[index][0]][2] + slicer.positions[slicer.triangles[index][1]][2] + slicer.positions[slicer.triangles[index][2]][2]) / 3;
	double value_of_manufacturing_require = (max_z - average_z) / (max_z - min_z);
	return value_of_manufacturing_require;
}

void HybridManufacturing::GetALLFragileVertex(SAMPLE_ON_BALL sampling)
{
	for (int ori = 0; ori < sampling.sample_points.size(); ori++)
	{
		// rotating the blocks and then slicing//
		vector<Eigen::MatrixXd> temp_V;
		temp_V.resize(V.rows());
		for (int i = 0; i < V.rows(); i++)
		{
			temp_V[i].resize(3, 1);
			temp_V[i](0, 0) = V.row(i).x();
			temp_V[i](1, 0) = V.row(i).y();
			temp_V[i](2, 0) = V.row(i).z();
		}
		Eigen::Vector3d vectorBefore(0, 0, 1);
		Eigen::Vector3d vectorAfter(sampling.sample_points[ori].x, sampling.sample_points[ori].y, sampling.sample_points[ori].z);
		Eigen::Matrix3d rotMatrix = Eigen::Quaterniond::FromTwoVectors(vectorBefore, vectorAfter).toRotationMatrix();
		for (int i = 0; i < V.rows(); i++)
			temp_V[i] = rotMatrix.inverse() * temp_V[i];

		vector<Vector3> normal_V;
		vector<bool> temp_fragile_V;
		normal_V.resize(V.rows());
		temp_fragile_V.resize(V.rows());
		for (int i = 0; i < V.rows(); i++)
		{
			normal_V[i].m_x = normal_V[i].m_y = normal_V[i].m_z = 0;
		}
		for (int i = 0; i < F.rows(); i++)
		{
			VEctor v1, v2, v3;
			v1[0] = temp_V[F(i, 0)](0, 0);
			v1[1] = temp_V[F(i, 0)](1, 0);
			v1[2] = temp_V[F(i, 0)](2, 0);
			v2[0] = temp_V[F(i, 1)](0, 0);
			v2[1] = temp_V[F(i, 1)](1, 0);
			v2[2] = temp_V[F(i, 1)](2, 0);
			v3[0] = temp_V[F(i, 2)](0, 0);
			v3[1] = temp_V[F(i, 2)](1, 0);
			v3[2] = temp_V[F(i, 2)](2, 0);
			Vector3 temp_normal;
			temp_normal = calculate_normal(v1, v2, v3);
			temp_normal.Normalized();
			normal_V[F(i, 0)] += temp_normal;
			normal_V[F(i, 1)] += temp_normal;
			normal_V[F(i, 2)] += temp_normal;
		}
		for (int i = 0; i < V.rows(); i++)
		{
			normal_V[i].Normalized();
			Vector3 temp_vectorBefore(0, 0, 1);
			if (Dot(normal_V[i], temp_vectorBefore) > 0.95)
				temp_fragile_V[i] = true;
			else
				temp_fragile_V[i] = false;
		}
		normal_all_ori_V.push_back(normal_V);
		is_fragile_V.push_back(temp_fragile_V);
	}
}

void HybridManufacturing::quick_sort(vector<TRiangle> &candidate_triangles, int left, int right, vector<int> &id_triangles, vector<double> &min_z_triangle, vector<VEctor> &min_z_point)
{
	if (left > right)
	{
		return;
	}

	double temp = min_z_triangle[left];
	TRiangle temp_2 = candidate_triangles[left];
	int temp_3 = id_triangles[left];
	VEctor temp_4 = min_z_point[left];

	int i = left;
	int j = right;
	int t;
	while (i != j)
	{

		while (i < j && min_z_triangle[j] >= temp)
		{
			j--;
		}

		while (i < j && min_z_triangle[i] <= temp)
		{
			i++;
		}

		if (i < j)
		{
			swap(candidate_triangles[i], candidate_triangles[j]);
			swap(id_triangles[i], id_triangles[j]);
			swap(min_z_triangle[i], min_z_triangle[j]);
			swap(min_z_point[i], min_z_point[j]);
		}
	}

	min_z_triangle[left] = min_z_triangle[i];
	min_z_triangle[i] = temp;
	candidate_triangles[left] = candidate_triangles[i];
	candidate_triangles[i] = temp_2;
	id_triangles[left] = id_triangles[i];
	id_triangles[i] = temp_3;
	min_z_point[left] = min_z_point[i];
	min_z_point[i] = temp_4;

	quick_sort(candidate_triangles, left, i - 1, id_triangles, min_z_triangle, min_z_point);
	quick_sort(candidate_triangles, i + 1, right, id_triangles, min_z_triangle, min_z_point);
}

void HybridManufacturing::Anticlockwise(vector<vector<int>> &real_cutting_plane_triangles, Slicer_2 all_slicer)
{

	for (int t = 0; t < real_cutting_plane_triangles.size(); t++)
	{
		double max_x = -MAX_D;
		int index_point;
		for (int i = 0; i < real_cutting_plane_triangles[t].size(); i++)
		{
			if (all_slicer.positions[real_cutting_plane_triangles[t][i]][0] > max_x)
			{
				max_x = all_slicer.positions[real_cutting_plane_triangles[t][i]][0];
				index_point = i;
			}
		}
		double d = (all_slicer.positions[real_cutting_plane_triangles[t][(index_point) % real_cutting_plane_triangles[t].size()]][0] - all_slicer.positions[real_cutting_plane_triangles[t][(index_point - 1 + real_cutting_plane_triangles[t].size()) % real_cutting_plane_triangles[t].size()]][0]) * (all_slicer.positions[real_cutting_plane_triangles[t][(index_point) % real_cutting_plane_triangles[t].size()]][1] - all_slicer.positions[real_cutting_plane_triangles[t][(index_point + 1) % real_cutting_plane_triangles[t].size()]][1]) - (all_slicer.positions[real_cutting_plane_triangles[t][(index_point) % real_cutting_plane_triangles[t].size()]][1] - all_slicer.positions[real_cutting_plane_triangles[t][(index_point - 1 + real_cutting_plane_triangles[t].size()) % real_cutting_plane_triangles[t].size()]][1]) * (all_slicer.positions[real_cutting_plane_triangles[t][(index_point) % real_cutting_plane_triangles[t].size()]][0] - all_slicer.positions[real_cutting_plane_triangles[t][(index_point + 1) % real_cutting_plane_triangles[t].size()]][0]);
		if (d > 0)
			reverse(real_cutting_plane_triangles[t].begin(), real_cutting_plane_triangles[t].end());
	}
}

int HybridManufacturing::isAnticlockwise(Point3f p1, Point3f p2, Point3f p3)
{
	//
	double ans = (p2.x - p1.x) * (p2.y - p3.y) - (p2.y - p1.y) * (p2.x - p3.x);
	// std::cout << ans << '\n';
	if (ans < 0) // is anticlockwise
		return 1;
	else if (ans == 0) // is line
		return 2;
	else if (ans > 0) // is clockwise
		return 3;
}

bool HybridManufacturing::checkHaveNoOtherPoint(Point3f p1, Point3f p2, Point3f p3, std::list<Point3f> pointList)
{
	std::list<Point3f>::iterator iter = pointList.begin();
	for (; iter != pointList.end(); iter++)
	{
		Point3f p = *iter;
		if (p == p1 || p == p2 || p == p3)
			continue;
		// if point p is in the triangle, (p,p1,p2),(p,p2,p3),(p,p3,p1) should be anticlockwise
		if (isAnticlockwise(p, p1, p2) == 1 && isAnticlockwise(p, p2, p3) == 1 && isAnticlockwise(p, p3, p1) == 1)
			return false;
	}
	return true;
}

void HybridManufacturing::detect_collision_with_printing_platform(int &index, vector<int> &candidate_nodes, vector<all_value> &all_calculated_value, vector<vector<cv::Point3d>> all_cut_layers, Eigen::Vector3d ori_now, nozzle the_nozzle)
{
	cv::Point3d center_point(0, 0, 0);
	for (int i = 0; i < V_bottom.size(); i++)
	{
		center_point.x += V_bottom[i].m_x;
		center_point.y += V_bottom[i].m_y;
		center_point.z += V_bottom[i].m_z;
	}
	center_point.x /= V_bottom.size();
	center_point.y /= V_bottom.size();
	center_point.z /= V_bottom.size();
	double circle_r = 25;
	if (ori_now.x() == 0 && ori_now.y() == 0 && ori_now.z() == 1)
		return;

	// 圆周上采样
	vector<Eigen::MatrixXd> base_points;
	int numPoints = 100;
	double angleStep = 2 * M_PI / double(numPoints); // 每个点的角度间隔

	for (int i = 0; i < numPoints; ++i)
	{
		double angle = i * angleStep;
		Eigen::MatrixXd point;
		point.resize(3, 1);
		point(0, 0) = center_point.x + circle_r * cos(angle); // 计算x坐标
		point(1, 0) = center_point.y + circle_r * sin(angle); // 计算y坐标
		point(2, 0) = center_point.z;
		base_points.push_back(point);
	}

	// 旋转
	Eigen::Vector3d vectorBefore(0, 0, 1);
	Eigen::Matrix3d rotMatrix = Eigen::Quaterniond::FromTwoVectors(vectorBefore, ori_now).toRotationMatrix();
	vector<vector<Eigen::MatrixXd>> rotate_all_cut_layers;
	rotate_all_cut_layers.resize(all_cut_layers.size());
	for (int i = 0; i < all_cut_layers.size(); i++)
	{
		rotate_all_cut_layers[i].resize(all_cut_layers[i].size());
		for (int j = 0; j < all_cut_layers[i].size(); j++)
		{
			rotate_all_cut_layers[i][j].resize(3, 1);
			rotate_all_cut_layers[i][j](0, 0) = all_cut_layers[i][j].x;
			rotate_all_cut_layers[i][j](1, 0) = all_cut_layers[i][j].y;
			rotate_all_cut_layers[i][j](2, 0) = all_cut_layers[i][j].z;
			// rotate_all_cut_layers[i][j] = rotMatrix.inverse() * rotate_all_cut_layers[i][j];
		}
	}
	for (int i = 0; i < base_points.size(); i++)
		base_points[i] = rotMatrix.inverse() * base_points[i];

	// 碰撞检测
	bool jud_collision = false;
	for (int i = 0; i < all_cut_layers.size(); i++)
	{
		cv::Point3d detect_point;
		for (int j = 0; j < all_cut_layers[i].size(); j++)
		{
			detect_point = cv::Point3d(rotate_all_cut_layers[i][j](0, 0), rotate_all_cut_layers[i][j](1, 0), rotate_all_cut_layers[i][j](2, 0));
			for (int k = 0; k < base_points.size(); k++)
			{
				if (base_points[k](2, 0) - detect_point.z <= 0)
				{
					jud_collision = false;
					continue;
				}
				if (base_points[k](2, 0) - detect_point.z > the_nozzle.nozzle__H_total)
				{
					jud_collision = true;
					break;
				}
				if (pow(base_points[k](0, 0) - detect_point.x, 2) + pow(base_points[k](1, 0) - detect_point.y, 2) - pow(the_nozzle.upper_surface_r, 2) < 0)
				{
					jud_collision = true;
					break;
				}
			}
			if (jud_collision == true)
				break;
		}
		if (jud_collision == true)
			break;
	}
	if (jud_collision == true)
	{
		candidate_nodes.erase(candidate_nodes.begin() + index);
		all_calculated_value.erase(all_calculated_value.begin() + index);
		index--;
		return;
	}
}

bool HybridManufacturing::Check_terminate(Slicer_2 slicer)
{
	double max_z;
	double min_z;
	for (int i = 0; i < slicer.triangles.size(); i++)
	{
		for (int j = 0; j < 3; j++)
		{
			if (i == 0 && j == 0)
			{
				max_z = min_z = slicer.positions[slicer.triangles[i][j]][2];
			}
			else
			{
				max_z = std::max(max_z, slicer.positions[slicer.triangles[i][j]][2]);
				min_z = std::min(min_z, slicer.positions[slicer.triangles[i][j]][2]);
			}
		}
	}
	if (max_z - min_z < 8)
		return true;
	else
		return false;
}

double HybridManufacturing::calculate_projected_area(Slicer_2 all_slicer, vector<vector<TRiangle>> all_furcation_of_blocks, vector<vector<cv::Point3d>> all_cut_layers)
{
	typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
	typedef K::Point_2 Point_2;
	typedef CGAL::Convex_hull_traits_adapter_2<K,
											   CGAL::Pointer_property_map<Point_2>::type>
		Convex_hull_traits_2;
	vector<vector<cv::Point3d>> vis_points;
	vis_points.resize(all_cut_layers.size());
	double avg_value = 0;

	for (int i = 0; i < all_cut_layers.size(); i++)
	{
		vector<Point_2> points_in_interface;
		points_in_interface.clear();
		double area_of_inter_face = 0;
		double area_of_projected_boundary = 0;
		vector<Point_2> points_in_furcation_of_blocks;
		points_in_furcation_of_blocks.clear();
		vector<Point_2> points_in_projected_boundary;
		points_in_projected_boundary.clear();

		for (int j = 0; j < all_cut_layers[i].size(); j++)
		{
			Point_2 temp_point(all_cut_layers[i][j].x, all_cut_layers[i][j].y);
			points_in_interface.push_back(temp_point);
		}
		for (int j = 0; j < points_in_interface.size(); j++)
		{
			area_of_inter_face += (points_in_interface[j].x() * points_in_interface[(j + 1) % points_in_interface.size()].y()) - (points_in_interface[(j + 1) % points_in_interface.size()].x() * points_in_interface[j].y());
		}
		area_of_inter_face /= 2;
		area_of_inter_face = abs(area_of_inter_face);

		for (int j = 0; j < all_furcation_of_blocks[i].size(); j++)
		{
			for (int k = 0; k < 3; k++)
			{
				Point_2 temp_point(all_slicer.positions[all_furcation_of_blocks[i][j][k]][0], all_slicer.positions[all_furcation_of_blocks[i][j][k]][1]);
				points_in_furcation_of_blocks.push_back(temp_point);
				// cv::Point3d temp_3d_point(temp_point.x(), temp_point.y(), 0);
				// vis_points[i].push_back(temp_3d_point);
			}
		}
		std::vector<std::size_t> indices(points_in_furcation_of_blocks.size()), out;
		std::iota(indices.begin(), indices.end(), 0);
		CGAL::convex_hull_2(indices.begin(), indices.end(), std::back_inserter(out),
							Convex_hull_traits_2(CGAL::make_property_map(points_in_furcation_of_blocks)));
		for (std::size_t t : out)
		{
			// std::cout << "points[" << t << "] = " << points_in_furcation_of_blocks[t] << std::endl;
			points_in_projected_boundary.push_back(points_in_furcation_of_blocks[t]);
			// cv::Point3d temp_3d_point(points_in_furcation_of_blocks[t].x(), points_in_furcation_of_blocks[t].y(), 0);
			// vis_points[i].push_back(temp_3d_point);
		}

		for (int j = 0; j < points_in_projected_boundary.size(); j++)
		{
			area_of_projected_boundary += (points_in_projected_boundary[j].x() * points_in_projected_boundary[(j + 1) % points_in_projected_boundary.size()].y()) - (points_in_projected_boundary[(j + 1) % points_in_projected_boundary.size()].x() * points_in_projected_boundary[j].y());
		}
		area_of_projected_boundary /= 2;
		area_of_projected_boundary = abs(area_of_projected_boundary);

		int n1 = points_in_projected_boundary.size(), n2 = points_in_interface.size();
		PolyIntersec polyint;
		MyPoint ps1[10000];
		MyPoint ps2[10000];
		for (int j = 0; j < n1; j++)
		{
			ps1[j].x = points_in_projected_boundary[j].x();
			ps1[j].y = points_in_projected_boundary[j].y();
		}
		for (int j = 0; j < n2; j++)
		{
			ps2[j].x = points_in_interface[j].x();
			ps2[j].y = points_in_interface[j].y();
		}
		double ans = polyint.intersectArea(ps1, n1, ps2, n2);
		ans = fabs(polyint.area(ps1, n1)) + fabs(polyint.area(ps2, n2)) - ans; // 容斥
		// cout << ans<<endl;

		avg_value += ((area_of_projected_boundary) / area_of_inter_face);
	}
	avg_value /= all_cut_layers.size();

	// Visual vis;
	// vis.generateModelForRendering_9(vis_points, ".\\vis\\projected");
	return avg_value;
}

all_value HybridManufacturing::GainMesh(Slicer_2 &slicer, vector<vector<cv::Point3d>> all_cut_layers, Eigen::Vector3d vector_after, int height_of_beam_search, int cont_number_of_queue, int index_of_pre_node, vector<int> all_cut_layers_dependency_layer, bool flag_is_continue_block, int id_continue)
{
	clock_t start_time, end_time;

	double sum_area = 0;
	double sum_value_of_manufacturing_require = 0;
	double sum_not_self_support_area = 0;
	// double area_threshold = 5.0;
	double ratio_of_area_threshold = 0.05; // 0.015 //0.15  //0.08   //0.03  //0.05
	double max_self_support_slope_angle;
	// cout << "#####" << height_of_beam_search << " " << id_continue << endl;
	// if(height_of_beam_search != 7 || id_continue != 1)s
	//	max_self_support_slope_angle = PI / 10;  //PI/3.6
	// else
	max_self_support_slope_angle = PI / 3.6; // PI/3.6
	double sum_area_threshold = 1000;		 // 200 //1000
	if (vector_after.x() == 0 && vector_after.y() == 0)
	{
		sum_area_threshold = 100000;
		ratio_of_area_threshold = 0.1;
	}

	/*if(height_of_beam_search == 4)
		ratio_of_area_threshold = 0.2;*/

	std::vector<bool> jud_triangle_have_been_added;
	double self_support_value;
	vector<vector<int>> cutting_plane_points;
	cutting_plane_points.resize(all_cut_layers.size());
	//////sort cut layers//////
	for (int i = 0; i < all_cut_layers.size(); i++)
		for (int j = i + 1; j < all_cut_layers.size(); j++)
		{

			if (all_cut_layers[i][0].z > all_cut_layers[j][0].z)
			{
				swap(all_cut_layers[i], all_cut_layers[j]);
				swap(all_cut_layers_dependency_layer[i], all_cut_layers_dependency_layer[j]);
			}
		}
	///////////////////////////

	// need rotate first//
	Eigen::Vector3d vectorBefore(0, 0, 1);
	Eigen::Matrix3d rotMatrix = Eigen::Quaterniond::FromTwoVectors(vectorBefore, vector_after).toRotationMatrix();
	for (int i = 0; i < slicer.positions.size(); i++)
	{
		Eigen::MatrixXd temp_V;
		temp_V.resize(3, 1);
		temp_V(0, 0) = slicer.positions[i][0];
		temp_V(1, 0) = slicer.positions[i][1];
		temp_V(2, 0) = slicer.positions[i][2];
		temp_V = rotMatrix.inverse() * temp_V;
		slicer.positions[i][0] = temp_V(0, 0);
		slicer.positions[i][1] = temp_V(1, 0);
		slicer.positions[i][2] = temp_V(2, 0);
	}
	/////////////////////

	//--------------------cut-----------------------//
	slicer.normal[0] = 0;
	slicer.normal[1] = 0;
	slicer.normal[2] = 1;
	// cout << endl << "pp" << endl;
	for (int i = 0; i < all_cut_layers.size(); i++)
	{
		slicer.origin[0] = 0;
		slicer.origin[1] = 0;
		slicer.origin[2] = all_cut_layers[i][0].z;
		slicer.cut();
	}
	// cout << endl << "qq" << endl;
	std::vector<TRiangle> ori_triangle = slicer.triangles;
	jud_triangle_have_been_added.resize(slicer.triangles.size());
	for (int i = 0; i < slicer.triangles.size(); i++)
		jud_triangle_have_been_added[i] = false;

	Slicer_2 all_slicer;

	all_slicer.positions = slicer.positions;
	all_slicer.triangles = ori_triangle;

	//--------------------save candidate_triangles-----------------------//
	int current_index = 0;
	std::vector<VEctor> min_z_point;
	std::vector<double> min_z_triangle;
	vector<int> index_of_min_point_in_triangle;
	std::vector<TRiangle> candidate_triangles;
	std::vector<int> id_triangles;
	std::vector<TRiangle> remove_triangles;
	candidate_triangles.clear();
	id_triangles.clear();
	double min_z_all_cut_layers = 9999999;
	for (int t = 0; t < all_cut_layers.size(); t++)
		min_z_all_cut_layers = min(min_z_all_cut_layers, all_cut_layers[t][0].z);
	for (int i = 0; i < all_slicer.triangles.size(); i++)
	{
		double temp_min_z_triangle = 9999999;
		for (int k = 0; k < 3; k++)
		{
			temp_min_z_triangle = min(all_slicer.positions[all_slicer.triangles[i][k]][2], temp_min_z_triangle);
		}
		if (temp_min_z_triangle + 0.001 >= min_z_all_cut_layers)
		{
			candidate_triangles.push_back(all_slicer.triangles[i]);
			id_triangles.push_back(i);
		}
	}
	for (int i = 0; i < candidate_triangles.size(); i++)
	{
		double temp_min_z_triangle = 9999999;
		VEctor temp_min_z_point = all_slicer.positions[candidate_triangles[i][0]];
		int temp_index_of_min_point_in_triangle;
		for (int k = 0; k < 3; k++)
		{
			if (all_slicer.positions[candidate_triangles[i][k]][2] < temp_min_z_triangle)
			{
				temp_min_z_triangle = all_slicer.positions[candidate_triangles[i][k]][2];
				temp_min_z_point = all_slicer.positions[candidate_triangles[i][k]];
				temp_index_of_min_point_in_triangle = candidate_triangles[i][k];
			}
		}
		min_z_triangle.push_back(temp_min_z_triangle);
		min_z_point.push_back(temp_min_z_point);
		index_of_min_point_in_triangle.push_back(temp_index_of_min_point_in_triangle);
	}

	quick_sort(candidate_triangles, 0, candidate_triangles.size() - 1, id_triangles, min_z_triangle, min_z_point);

	////--------------------save OPP_triangles one by one-----------------------//
	vector<vector<TRiangle>> all_furcation_of_blocks;
	all_furcation_of_blocks.resize(all_cut_layers.size());
	vector<int> id_of_furcation_of_blocks;
	std::vector<double> boundary_of_top_layer;
	vector<int> save_current_index;
	for (int t = 0; t < all_cut_layers.size(); t++)
	{
		current_index = 0; // 该方式也许比较慢
		double boundary_bottom = 999999, boundary_left = 999999, boundary_top = -999999, boundary_right = -999999;
		for (int i = 0; i < all_cut_layers[t].size(); i++)
		{
			boundary_top = std::max(boundary_top, all_cut_layers[t][i].y);
			boundary_bottom = std::min(boundary_bottom, all_cut_layers[t][i].y);
			boundary_right = std::max(boundary_right, all_cut_layers[t][i].x);
			boundary_left = std::min(boundary_left, all_cut_layers[t][i].x);
		}
		cv::Point2d current_triangle_point;
		cv::Point2d current_layer_point;
		for (; current_index < candidate_triangles.size(); current_index++)
		{
			if (abs(min_z_triangle[current_index] - all_cut_layers[t][0].z) > 0.001)
			{
				if (min_z_triangle[current_index] > all_cut_layers[t][0].z)
					break;
			}
			else
			{
				bool jud_is_boundary_point = false;
				// for (int i = 0; i < all_cut_layers[t].size(); i++) {
				int cont_inside_boundary = 0;
				for (int k = 0; k < 3; k++)
				{
					if (all_slicer.positions[candidate_triangles[current_index][k]][0] + 0.01 >= boundary_left && all_slicer.positions[candidate_triangles[current_index][k]][0] - 0.01 <= boundary_right && all_slicer.positions[candidate_triangles[current_index][k]][1] + 0.01 >= boundary_bottom && all_slicer.positions[candidate_triangles[current_index][k]][1] - 0.01 <= boundary_top)
					{
						cont_inside_boundary++;
					}
				}
				if (cont_inside_boundary >= 2)
					jud_is_boundary_point = true;
				/*if (min_z_point[current_index][0] + 0.01 >= boundary_left && min_z_point[current_index][0] - 0.01 <= boundary_right &&
					min_z_point[current_index][1] + 0.01 >= boundary_bottom && min_z_point[current_index][1] - 0.01 <= boundary_top) {
					jud_is_boundary_point = true;
					break;
				}*/
				//}
				if (jud_is_boundary_point == true)
				{
					for (int j = 0; j < all_cut_layers[t].size(); j++)
					{
						current_layer_point.x = all_cut_layers[t][j].x;
						current_layer_point.y = all_cut_layers[t][j].y;
						current_triangle_point.x = min_z_point[current_index][0];
						current_triangle_point.y = min_z_point[current_index][1];
						if (distance_2d(current_layer_point, current_triangle_point) < 0.2)
						{ // 4.0
							// if (height_of_beam_search != 2)
							remove_triangles.push_back(candidate_triangles[current_index]);
							id_of_furcation_of_blocks.push_back(t);
							all_furcation_of_blocks[t].push_back(candidate_triangles[current_index]);

							Vector3 face_normal = calculate_normal(all_slicer.positions[remove_triangles[remove_triangles.size() - 1][0]], all_slicer.positions[remove_triangles[remove_triangles.size() - 1][1]], all_slicer.positions[remove_triangles[remove_triangles.size() - 1][2]]);
							face_normal.Normalized();
							Vector3 base_normal(0, 0, 1);
							int jud_self_support = (Dot(face_normal, base_normal) + sin(max_self_support_slope_angle) >= 0);
							double current_area = calculate_area(all_slicer.positions[remove_triangles[remove_triangles.size() - 1][0]], all_slicer.positions[remove_triangles[remove_triangles.size() - 1][1]], all_slicer.positions[remove_triangles[remove_triangles.size() - 1][2]]);
							sum_area += current_area;
							// sum_value_of_manufacturing_require += current_area * calculate_manufacturing_require(all_slicer, remove_triangles);
							if (jud_self_support == false)
							{
								sum_not_self_support_area += current_area;
								if (sum_not_self_support_area / sum_area > ratio_of_area_threshold && sum_area >= sum_area_threshold)
								{
									all_value temp_all_value;
									temp_all_value.value_of_self_support = 0;
									return temp_all_value;
								}
							}

							save_current_index.push_back(current_index);
							jud_triangle_have_been_added[id_triangles[current_index]] = true;
							break;
						}
					}
				}
			}
		}
	}
	current_index = 0; // current_index = 0;
	// std::cout << "**" << current_index << " " << candidate_triangles.size() << " " << remove_triangles.size() << endl;

	start_time = clock();
	int index_current_cut_layer = 1;
	while (1)
	{
		bool flag_break = false;
		for (int i = 0; i < remove_triangles.size(); i++)
		{
			if (jud_triangle_have_been_added[id_triangles[current_index]] == false && min_z_triangle[current_index] >= min_z_triangle[save_current_index[i]])
			{
				for (int j = 0; j < 3; j++)
				{
					for (int k = 0; k < 3; k++)
					{
						if (candidate_triangles[current_index][j] == remove_triangles[i][k])
						{
							// if(height_of_beam_search != 2)
							remove_triangles.push_back(candidate_triangles[current_index]);
							id_of_furcation_of_blocks.push_back(id_of_furcation_of_blocks[i]);
							all_furcation_of_blocks[id_of_furcation_of_blocks[i]].push_back(candidate_triangles[current_index]);
							//////all_furcation_of_blocks[t].push_back(candidate_triangles[current_index]);

							Vector3 face_normal = calculate_normal(all_slicer.positions[remove_triangles[remove_triangles.size() - 1][0]], all_slicer.positions[remove_triangles[remove_triangles.size() - 1][1]], all_slicer.positions[remove_triangles[remove_triangles.size() - 1][2]]);
							face_normal.Normalized();
							Vector3 base_normal(0, 0, 1);
							int jud_self_support = (Dot(face_normal, base_normal) + sin(max_self_support_slope_angle) >= 0);
							double current_area = calculate_area(all_slicer.positions[remove_triangles[remove_triangles.size() - 1][0]], all_slicer.positions[remove_triangles[remove_triangles.size() - 1][1]], all_slicer.positions[remove_triangles[remove_triangles.size() - 1][2]]);
							sum_area += current_area;

							if (jud_self_support == false)
							{
								sum_not_self_support_area += current_area;
								if (sum_not_self_support_area / sum_area > ratio_of_area_threshold && sum_area >= sum_area_threshold)
								{
									all_value temp_all_value;
									temp_all_value.value_of_self_support = 0;
									return temp_all_value;
								}
							}

							save_current_index.push_back(current_index);
							jud_triangle_have_been_added[id_triangles[current_index]] = true;
							flag_break = true;
							break;
						}
					}
					if (flag_break == true)
						break;
				}
			}
			if (flag_break == true)
				break;
		}
		current_index++;
		if (current_index >= candidate_triangles.size())
			break;
	}
	end_time = clock();
	// add remaining face, set a distance threshold Dis, only a face less Dis from other cut layers and no other dependency layer exist, the face is consider to remaining face
	/*double Dis = dh * 4;
	for (int i = 0; i < jud_triangle_have_been_added.size(); i++) {
		if (jud_triangle_have_been_added[i] == false) {
			int id_layer;
			double min_dis = 99999999;
			for (int t = 0; t < all_cut_layers.size(); t++) {
				if (all_cut_layers[t][0].z - slicer.positions[slicer.triangles[i][0]][2] <= 2 * dh && all_cut_layers[t][0].z - slicer.positions[slicer.triangles[i][0]][2] >= -0.001) {
					for (int j = 0; j < all_cut_layers[t].size(); j++) {
						cv::Point3d current_triangle_point(slicer.positions[slicer.triangles[i][0]][0], slicer.positions[slicer.triangles[i][0]][1], slicer.positions[slicer.triangles[i][0]][2]);
						cv::Point3d current_layer_point(all_cut_layers[t][j].x, all_cut_layers[t][j].y, all_cut_layers[t][0].z);

						double distance = distance_3d(current_triangle_point, current_layer_point);
						if (distance < min_dis) {
							min_dis = distance;
							id_layer = t;
						}
					}
				}
			}
			if (min_dis < Dis && all_cut_layers_dependency_layer[id_layer] == 0)
				remove_triangles.push_back(slicer.triangles[i]);
		}
	}*/

	all_value all_calculated_value;
	all_calculated_value.number_of_remaining_face = all_slicer.triangles.size() - remove_triangles.size();
	all_slicer.triangles = remove_triangles;

	all_calculated_value.value_of_projected = calculate_projected_area(all_slicer, all_furcation_of_blocks, all_cut_layers);

	double min_z = 9999999, max_z = -9999999;
	double value_of_manufacturing_require;
	for (int i = 0; i < all_slicer.positions.size(); i++)
	{
		min_z = min(min_z, all_slicer.positions[i][2]);
		max_z = max(max_z, all_slicer.positions[i][2]);
	}
	for (int i = 0; i < all_slicer.triangles.size(); i++)
	{
		double current_area = calculate_area(all_slicer.positions[all_slicer.triangles[i][0]], all_slicer.positions[all_slicer.triangles[i][1]], all_slicer.positions[all_slicer.triangles[i][2]]);
		value_of_manufacturing_require = calculate_manufacturing_require(all_slicer, i, min_z, max_z);
		sum_value_of_manufacturing_require += value_of_manufacturing_require * current_area;
	}

	all_calculated_value.value_of_self_support = 1 - sum_not_self_support_area / sum_area;
	all_calculated_value.large_base = sum_value_of_manufacturing_require / sum_area;

	// cout << "()()()(" << double(end_time - start_time) / CLOCKS_PER_SEC << endl;
	return all_calculated_value;
}

void HybridManufacturing::CutMesh(vector<vector<cv::Point3d>> all_layers, vector<vector<cv::Point3d>> all_layers_contain, vector<vector<cv::Point3d>> all_cut_layers, Eigen::Vector3d vector_after, int height_of_beam_search, int cont_number_of_queue, int index_of_pre_node, vector<int> all_cut_layers_dependency_layer, bool &jud_outer_beam_search_terminate, vector<TRiangle> &current_remove_triangles, Slicer_2 &current_slicer, bool judge_continue_additive, bool flag_is_continue_block, int pre_cont_number_of_queue, vector<bool> &jud_error, int id_node, int id_continue, vector<int> flag_cut_layers_is_hole)
{
	bool using_solid_model = true;

	std::vector<bool> jud_triangle_have_been_added;
	vector<vector<int>> cutting_plane_points;
	cutting_plane_points.resize(all_cut_layers.size());
	vector<vector<pair<int, int>>> cutting_plane_edges;
	cutting_plane_edges.resize(all_cut_layers.size());
	//////sort cut layers//////
	map<int, int> follow_index;
	follow_index.clear();
	for (int i = 0; i < all_cut_layers.size(); i++)
		follow_index.insert({i, i});
	for (int i = 0; i < all_cut_layers.size(); i++)
		for (int j = i + 1; j < all_cut_layers.size(); j++)
		{
			if (all_cut_layers[i][0].z > all_cut_layers[j][0].z)
			{
				swap(all_cut_layers[i], all_cut_layers[j]);
				swap(follow_index[i], follow_index[j]);
				swap(flag_cut_layers_is_hole[i], flag_cut_layers_is_hole[j]);
				swap(all_cut_layers_dependency_layer[i], all_cut_layers_dependency_layer[j]);
			}
		}
	///////////////////////////
	// cout << "a" << id_continue;
	Slicer_2 slicer;
	if (flag_is_continue_block == false)
	{
		/*vector<std::array<double, 3>> V3;
		vector<std::array<int, 3>> F3;
		vector<std::array<double, 3>> N3;
		ifstream ifs(file_name + "-" + to_string(height_of_beam_search - 1) + "_" + to_string(index_of_pre_node) + ".stl");
		igl::read_stl_ascii(ifs, V3, F3, N3);
		Eigen::MatrixXd V4(V3.size(), 3);
		Eigen::MatrixXi F4(F3.size(), 3);
		Eigen::MatrixXd N4(N3.size(), 3);
		for (int i = 0; i < V3.size(); i++)
			for (int j = 0; j < 3; j++)
				V4(i, j) = V3[i][j];
		for (int i = 0; i < F3.size(); i++)
			for (int j = 0; j < 3; j++)
				F4(i, j) = F3[i][j];
		cout << V3.size() << "()()()()()" << endl;
		igl::writeOBJ(file_name + "-" + to_string(height_of_beam_search - 1) + "_" + to_string(index_of_pre_node) + ".obj", V4, F4);*/

		//  create a geometry tesselation object
		slicer.load(file_name + "-" + to_string(height_of_beam_search - 1) + "_" + to_string(index_of_pre_node) + ".obj");
		cout << "&" << endl;
	}
	else
	{
		/*Eigen::MatrixXd V3;
		Eigen::MatrixXi F3;
		Eigen::MatrixXi N3;
		FILE* fp;
		fp = fopen((file_name + "-" + to_string(height_of_beam_search - 1) + "_" + to_string(index_of_pre_node) + "_" + to_string(id_continue - 1) + "_subblock.stl").c_str(), "r");*/
		// igl::readSTL(fp, V3, F3, N3);
		// igl::writeOBJ(file_name + "-" + to_string(height_of_beam_search - 1) + "_" + to_string(index_of_pre_node) + "_" + to_string(id_continue - 1) + "_subblock.obj", V3, F3);
		slicer.load(file_name + "-" + to_string(height_of_beam_search - 1) + "_" + to_string(index_of_pre_node) + "_" + to_string(id_continue - 1) + "_subblock.obj");
		cout << "*" << endl;
	}

	if (flag_is_continue_block == true)
	{
		height_of_beam_search--;
		cont_number_of_queue = pre_cont_number_of_queue;
	}

	// need rotate first//
	Eigen::Vector3d vectorBefore(0, 0, 1);
	Eigen::Matrix3d rotMatrix = Eigen::Quaterniond::FromTwoVectors(vectorBefore, vector_after).toRotationMatrix();
	for (int i = 0; i < slicer.positions.size(); i++)
	{
		Eigen::MatrixXd temp_V;
		temp_V.resize(3, 1);
		temp_V(0, 0) = slicer.positions[i][0];
		temp_V(1, 0) = slicer.positions[i][1];
		temp_V(2, 0) = slicer.positions[i][2];
		temp_V = rotMatrix.inverse() * temp_V;
		slicer.positions[i][0] = temp_V(0, 0);
		slicer.positions[i][1] = temp_V(1, 0);
		slicer.positions[i][2] = temp_V(2, 0);
	}
	/////////////////////

	// layer visualization//
	rotMatrix = Eigen::Quaterniond::FromTwoVectors(vector_after, vectorBefore).toRotationMatrix();
	// if (height_of_beam_search != 2)
	for (int i = 0; i < all_layers.size(); i++)
	{
		for (int j = 0; j < all_layers[i].size(); j++)
		{
			Eigen::MatrixXd temp_V;
			temp_V.resize(3, 1);
			temp_V(0, 0) = all_layers[i][j].x;
			temp_V(1, 0) = all_layers[i][j].y;
			temp_V(2, 0) = all_layers[i][j].z;
			temp_V = rotMatrix.inverse() * temp_V;
			all_layers[i][j].x = temp_V(0, 0);
			all_layers[i][j].y = temp_V(1, 0);
			all_layers[i][j].z = temp_V(2, 0);
		}
		for (int j = 0; j < all_layers_contain[i].size(); j++)
		{
			Eigen::MatrixXd temp_V;
			temp_V.resize(3, 1);
			temp_V(0, 0) = all_layers_contain[i][j].x;
			temp_V(1, 0) = all_layers_contain[i][j].y;
			temp_V(2, 0) = all_layers_contain[i][j].z;
			temp_V = rotMatrix.inverse() * temp_V;
			all_layers_contain[i][j].x = temp_V(0, 0);
			all_layers_contain[i][j].y = temp_V(1, 0);
			all_layers_contain[i][j].z = temp_V(2, 0);
		}
	}
	cout << "b" << endl;
	Visual Vis;
	cout << id_continue << endl;
	// Vis.generateModelForRendering_5(all_layers, all_layers_contain, height_of_beam_search, cont_number_of_queue, file_name, index_of_pre_node,judge_continue_additive, id_continue);
	string vis_file(file_name + "-" + to_string(height_of_beam_search) + "_" + to_string(cont_number_of_queue) + "(" + to_string(index_of_pre_node) + ")_Layer");
	string vis_file_contain(file_name + "-" + to_string(height_of_beam_search) + "_" + to_string(cont_number_of_queue) + "(" + to_string(index_of_pre_node) + ")_Layer_contain");
	Point3ff current_orientation(vector_after.x(), vector_after.y(), vector_after.z());
	// visualize_layers(vis_file, current_orientation, judge_continue_additive, id_continue);
	if (open_vis_stair_case == true)
		visualize_layers_stair_case(vis_file, vis_file_contain, current_orientation, judge_continue_additive, id_continue);
	////////////////////////
	cout << "B" << endl;
	//--------------------cut-----------------------//
	slicer.normal[0] = 0;
	slicer.normal[1] = 0;
	slicer.normal[2] = 1;

	slicer.jud_plane.resize(slicer.triangles.size());
	for (int i = 0; i < slicer.triangles.size(); i++)
		slicer.jud_plane[i] = false;
	clock_t start_time_3, end_time_3;

	for (int i = 0; i < all_cut_layers.size(); i++)
	{
		slicer.origin[0] = 0;
		slicer.origin[1] = 0;
		slicer.origin[2] = all_cut_layers[i][0].z;
		slicer.cut();
		// cout << "ok";
	}

	std::vector<TRiangle> ori_triangle = slicer.triangles;
	jud_triangle_have_been_added.resize(slicer.triangles.size());
	for (int i = 0; i < slicer.triangles.size(); i++)
		jud_triangle_have_been_added[i] = false;

	Slicer_2 all_slicer;

	all_slicer.positions = slicer.positions;
	all_slicer.triangles = ori_triangle;
	//--------------------save candidate_triangles-----------------------//
	int current_index = 0;
	std::vector<VEctor> min_z_point;
	std::vector<double> min_z_triangle;
	vector<int> index_of_min_point_in_triangle;
	std::vector<TRiangle> candidate_triangles;
	std::vector<int> id_candidate_triangles;
	std::vector<int> id_triangles;
	std::vector<TRiangle> remove_triangles;
	std::vector<int> id_remove_triangles;
	candidate_triangles.clear();
	id_triangles.clear();
	double min_z_all_cut_layers = 9999999;
	for (int t = 0; t < all_cut_layers.size(); t++)
		min_z_all_cut_layers = min(min_z_all_cut_layers, all_cut_layers[t][0].z);
	start_time_3 = clock();
	for (int i = 0; i < all_slicer.triangles.size(); i++)
	{
		double temp_min_z_triangle = 9999999;
		for (int k = 0; k < 3; k++)
		{
			temp_min_z_triangle = min(all_slicer.positions[all_slicer.triangles[i][k]][2], temp_min_z_triangle);
		}
		if (temp_min_z_triangle + 0.001 >= min_z_all_cut_layers)
		{
			candidate_triangles.push_back(all_slicer.triangles[i]);
			id_candidate_triangles.push_back(i);
			for (int k = 0; k < slicer.triangles.size(); k++)
				if (slicer.triangles[k] == all_slicer.triangles[i])
				{
					id_triangles.push_back(k);
					break;
				}
		}
	}
	end_time_3 = clock();
	for (int i = 0; i < candidate_triangles.size(); i++)
	{
		double temp_min_z_triangle = 9999999;
		VEctor temp_min_z_point = all_slicer.positions[candidate_triangles[i][0]];
		int temp_index_of_min_point_in_triangle;
		for (int k = 0; k < 3; k++)
		{
			if (all_slicer.positions[candidate_triangles[i][k]][2] < temp_min_z_triangle)
			{
				temp_min_z_triangle = all_slicer.positions[candidate_triangles[i][k]][2];
				temp_min_z_point = all_slicer.positions[candidate_triangles[i][k]];
				temp_index_of_min_point_in_triangle = k;
			}
		}
		min_z_triangle.push_back(temp_min_z_triangle);
		min_z_point.push_back(temp_min_z_point);
		index_of_min_point_in_triangle.push_back(temp_index_of_min_point_in_triangle);
	}

	cout << "()()()(" << double(end_time_3 - start_time_3) / CLOCKS_PER_SEC << endl;
	clock_t start_time_2, end_time_2;
	start_time_2 = clock();
	for (int i = 0; i < candidate_triangles.size(); i++)
	{
		for (int j = i + 1; j < candidate_triangles.size(); j++)
		{
			if (min_z_triangle[i] > min_z_triangle[j])
			{
				swap(candidate_triangles[i], candidate_triangles[j]);
				swap(id_candidate_triangles[i], id_candidate_triangles[j]);
				swap(id_triangles[i], id_triangles[j]);
				swap(min_z_triangle[i], min_z_triangle[j]);
				swap(min_z_point[i], min_z_point[j]);
				swap(index_of_min_point_in_triangle[i], index_of_min_point_in_triangle[j]);
			}
		}
	}
	end_time_2 = clock();
	cout << "()()()(" << double(end_time_2 - start_time_2) / CLOCKS_PER_SEC << endl;
	// slicer.triangles = candidate_triangles;

	//--------------------save OPP_triangles one by one-----------------------//
	std::vector<double> boundary_of_top_layer;
	int num_boundary_remove_faces;
	vector<int> save_current_index;
	for (int t = 0; t < all_cut_layers.size(); t++)
	{
		current_index = 0; // 该方式也许比较慢
		double boundary_bottom = 999999, boundary_left = 999999, boundary_top = -999999, boundary_right = -999999;
		Point_2 *points = new Point_2[all_cut_layers[t].size()];
		for (int i = 0; i < all_cut_layers[t].size(); i++)
		{
			boundary_top = std::max(boundary_top, all_cut_layers[t][i].y);
			boundary_bottom = std::min(boundary_bottom, all_cut_layers[t][i].y);
			boundary_right = std::max(boundary_right, all_cut_layers[t][i].x);
			boundary_left = std::min(boundary_left, all_cut_layers[t][i].x);
			Point_2 temp_point(all_cut_layers[t][i].x, all_cut_layers[t][i].y);
			points[i] = temp_point;
		}
		cv::Point2d current_triangle_point;
		cv::Point2d current_layer_point;
		for (; current_index < candidate_triangles.size(); current_index++)
		{
			if (abs(min_z_triangle[current_index] - all_cut_layers[t][0].z) > 0.0001)
			{
				if (min_z_triangle[current_index] > all_cut_layers[t][0].z)
					break;
			}
			else
			{
				bool jud_is_boundary_point = true;
				// for (int i = 0; i < all_cut_layers[t].size(); i++) {
				/*int cont_inside_boundary = 0;
				for (int k = 0; k < 3; k++) {
					if (all_slicer.positions[candidate_triangles[current_index][k]][0] + 0.01 >= boundary_left && all_slicer.positions[candidate_triangles[current_index][k]][0] - 0.01 <= boundary_right
						&& all_slicer.positions[candidate_triangles[current_index][k]][1] + 0.01 >= boundary_bottom && all_slicer.positions[candidate_triangles[current_index][k]][1] - 0.01 <= boundary_top) {
						cont_inside_boundary++;
					}
				}
				if (cont_inside_boundary >= 2)
					jud_is_boundary_point = true;*/
				/*if (min_z_point[current_index][0] + 0.01 >= boundary_left && min_z_point[current_index][0] - 0.01 <= boundary_right &&
					min_z_point[current_index][1] + 0.01 >= boundary_bottom && min_z_point[current_index][1] - 0.01 <= boundary_top) {
					jud_is_boundary_point = true;
					break;
				}*/
				//}

				if (jud_is_boundary_point == true)
				{
					current_triangle_point = cv::Point2d(min_z_point[current_index][0], min_z_point[current_index][1]);
					if (check_inside_2(Point_2(current_triangle_point.x, current_triangle_point.y), points, points + all_cut_layers[t].size(), K()))
					{
						int cont_cutting_points = 0;
						int temp_left = 0, temp_right = 0;
						for (int k = 0; k < 3; k++)
						{
							if (abs(all_cut_layers[t][0].z - all_slicer.positions[candidate_triangles[current_index][k]][2]) < 0.0001)
							{
								cont_cutting_points++;
								if (cont_cutting_points == 1)
									temp_left = candidate_triangles[current_index][k];
								else if (cont_cutting_points == 2)
								{
									temp_right = candidate_triangles[current_index][k];
									cutting_plane_edges[t].push_back(make_pair(temp_left, temp_right));
									break;
								}
								// cutting_plane_points[t].push_back(candidate_triangles[current_index][index_of_min_point_in_triangle[current_index]]);
							}
						}
						remove_triangles.push_back(candidate_triangles[current_index]);
						id_remove_triangles.push_back(id_candidate_triangles[current_index]);
						/*if (height_of_beam_search == 2 && id_triangles[current_index] == 10297)
							cout << "********************************** " << all_cut_layers[t][0].z << endl;*/

						save_current_index.push_back(current_index);
						jud_triangle_have_been_added[id_triangles[current_index]] = true;
					}
					else
					{
						for (int j = 0; j < all_cut_layers[t].size(); j++)
						{
							current_layer_point.x = all_cut_layers[t][j].x;
							current_layer_point.y = all_cut_layers[t][j].y;
							current_triangle_point.x = min_z_point[current_index][0];
							current_triangle_point.y = min_z_point[current_index][1];
							// if ((jud_is_boundary_point == false && distance_2d(current_layer_point, current_triangle_point) < 0.002) || (jud_is_boundary_point == true)) {   //0.002
							if (distance_2d(current_layer_point, current_triangle_point) < 0.1)
							{ // 4.0
								int cont_cutting_points = 0;
								int temp_left = 0, temp_right = 0;
								for (int k = 0; k < 3; k++)
								{
									if (abs(all_cut_layers[t][0].z - all_slicer.positions[candidate_triangles[current_index][k]][2]) < 0.0001)
									{
										cont_cutting_points++;
										if (cont_cutting_points == 1)
											temp_left = candidate_triangles[current_index][k];
										else if (cont_cutting_points == 2)
										{
											temp_right = candidate_triangles[current_index][k];
											cutting_plane_edges[t].push_back(make_pair(temp_left, temp_right));
											break;
										}
										// cutting_plane_points[t].push_back(candidate_triangles[current_index][index_of_min_point_in_triangle[current_index]]);
									}
								}
								remove_triangles.push_back(candidate_triangles[current_index]);
								id_remove_triangles.push_back(id_candidate_triangles[current_index]);
								/*if (height_of_beam_search == 2 && id_triangles[current_index] == 10297)
									cout << "********************************** " << all_cut_layers[t][0].z << endl;*/

								save_current_index.push_back(current_index);
								jud_triangle_have_been_added[id_triangles[current_index]] = true;
								break;
							}
						}
					}
				}
			}
		}
	}
	current_index = 0;
	num_boundary_remove_faces = remove_triangles.size();

	int index_current_cut_layer = 1;
	while (1)
	{
		bool flag_break = false;
		for (int i = 0; i < remove_triangles.size(); i++)
		{
			for (int j = 0; j < 3; j++)
			{
				for (int k = 0; k < 3; k++)
				{
					if (jud_triangle_have_been_added[id_triangles[current_index]] == false && candidate_triangles[current_index][j] == remove_triangles[i][k] && min_z_triangle[current_index] >= min_z_triangle[save_current_index[i]])
					{
						// if (height_of_beam_search != 2)
						remove_triangles.push_back(candidate_triangles[current_index]);
						id_remove_triangles.push_back(id_candidate_triangles[current_index]);
						save_current_index.push_back(current_index);
						jud_triangle_have_been_added[id_triangles[current_index]] = true;
						flag_break = true;
						break;
					}
				}
				if (flag_break == true)
					break;
			}
			if (flag_break == true)
				break;
		}
		current_index++;
		if (current_index >= candidate_triangles.size())
			break;
	}

	// add remaining face, set a distance threshold Dis, only a face less Dis from other cut layers and no other dependency layer exist, the face is consider to remaining face
	double Dis = dh * 2;
	for (int i = 0; i < jud_triangle_have_been_added.size(); i++)
	{
		if (jud_triangle_have_been_added[i] == false)
		{
			int id_layer;
			double min_dis = 99999999;
			for (int t = 0; t < all_cut_layers.size(); t++)
			{
				if (all_cut_layers[t][0].z - slicer.positions[slicer.triangles[i][0]][2] <= 2 * dh && all_cut_layers[t][0].z - slicer.positions[slicer.triangles[i][0]][2] >= -0.001)
				{
					for (int j = 0; j < all_cut_layers[t].size(); j++)
					{
						cv::Point3d current_triangle_point(slicer.positions[slicer.triangles[i][0]][0], slicer.positions[slicer.triangles[i][0]][1], slicer.positions[slicer.triangles[i][0]][2]);
						cv::Point3d current_layer_point(all_cut_layers[t][j].x, all_cut_layers[t][j].y, all_cut_layers[t][0].z);

						double distance = distance_3d(current_triangle_point, current_layer_point);
						if (distance < min_dis)
						{
							min_dis = distance;
							id_layer = t;
						}
					}
				}
			}
			if (min_dis < Dis && all_cut_layers_dependency_layer[id_layer] == 0)
			{
				remove_triangles.push_back(slicer.triangles[i]);
				id_remove_triangles.push_back(i);
				// jud_triangle_have_been_added[i] = true;
			}
		}
	}

	/////////////////////////删除残余面片////////////////////////////
	// 建立面片邻接关系
	vector<vector<int>> adjacent_faces(slicer.triangles.size());
	for (int i = 0; i < slicer.triangles.size(); i++)
	{
		if (jud_triangle_have_been_added[i] == false)
		{
			for (int j = i + 1; j < slicer.triangles.size(); j++)
			{
				if (jud_triangle_have_been_added[j] == false)
				{
					int cont_same_point = 0;
					for (int k = 0; k < 3; k++)
					{
						for (int l = 0; l < 3; l++)
						{
							if (slicer.triangles[i][k] == slicer.triangles[j][l])
							{
								cont_same_point++;
								break;
							}
						}
					}
					if (cont_same_point >= 2)
					{
						adjacent_faces[i].push_back(j);
						adjacent_faces[j].push_back(i);
					}
				}
			}
		}
	}
	// 分区
	bool *visited = new bool[slicer.triangles.size()];
	for (int i = 0; i < slicer.triangles.size(); i++)
		visited[i] = false;
	vector<vector<int>> q_faces(0);
	vector<vector<int>> save_faces(0);
	for (int i = 0; i < slicer.triangles.size(); i++)
	{
		if (jud_triangle_have_been_added[i] == false && visited[i] == false)
		{
			visited[i] = true;
			vector<int> q_temp;
			q_faces.push_back(q_temp);
			save_faces.push_back(q_temp);
			q_faces[q_faces.size() - 1].push_back(i);
			save_faces[save_faces.size() - 1].push_back(i);
			while (!q_faces[q_faces.size() - 1].empty())
			{
				int u_faces = q_faces[q_faces.size() - 1][q_faces[q_faces.size() - 1].size() - 1];
				q_faces[q_faces.size() - 1].pop_back();
				for (int j = 0; j < adjacent_faces[u_faces].size(); j++)
				{
					if (!visited[adjacent_faces[u_faces][j]] && jud_triangle_have_been_added[adjacent_faces[u_faces][j]] == false)
					{
						visited[adjacent_faces[u_faces][j]] = true;
						q_faces[q_faces.size() - 1].push_back(adjacent_faces[u_faces][j]);
						save_faces[save_faces.size() - 1].push_back(adjacent_faces[u_faces][j]);
					}
				}
			}
		}
	}
	// 移除残余面
	for (int i = 0; i < save_faces.size(); i++)
	{
		cout << "&&*&& " << save_faces[i].size() << endl;
		if (save_faces[i].size() < 80)
		{
			for (int j = 0; j < save_faces[i].size(); j++)
			{
				remove_triangles.push_back(slicer.triangles[save_faces[i][j]]);
				id_remove_triangles.push_back(save_faces[i][j]);
				jud_triangle_have_been_added[save_faces[i][j]] = true;
			}
		}
	}
	//////////////////////////////////////////////////////////////////

	cout << "c" << endl;
	// current_remove_triangles = remove_triangles;
	if (height_of_beam_search != 1)
	{
		string str_contact_faces;
		if (flag_is_continue_block == false)
			str_contact_faces = file_name + "-" + to_string(height_of_beam_search - 1) + "_" + to_string(index_of_pre_node) + "_contact_faces.txt";
		else
			str_contact_faces = file_name + "-" + to_string(height_of_beam_search) + "_" + to_string(index_of_pre_node) + "_" + to_string(id_continue - 1) + "_subblock_contact_faces.txt";
		ifstream load_contact_faces(str_contact_faces);
		int num_contact_faces;
		vector<int> id_last_contact_faces;
		vector<Vector3> last_contact_faces_normals;
		load_contact_faces >> num_contact_faces;
		id_last_contact_faces.resize(num_contact_faces);
		for (int j = 0; j < num_contact_faces; j++)
		{
			load_contact_faces >> id_last_contact_faces[j];
			Vector3 v1, v2, v3;
			v1.m_x = slicer.positions[slicer.triangles[id_last_contact_faces[j]][0]][0];
			v1.m_y = slicer.positions[slicer.triangles[id_last_contact_faces[j]][0]][1];
			v1.m_z = slicer.positions[slicer.triangles[id_last_contact_faces[j]][0]][2];
			v2.m_x = slicer.positions[slicer.triangles[id_last_contact_faces[j]][1]][0];
			v2.m_y = slicer.positions[slicer.triangles[id_last_contact_faces[j]][1]][1];
			v2.m_z = slicer.positions[slicer.triangles[id_last_contact_faces[j]][1]][2];
			v3.m_x = slicer.positions[slicer.triangles[id_last_contact_faces[j]][2]][0];
			v3.m_y = slicer.positions[slicer.triangles[id_last_contact_faces[j]][2]][1];
			v3.m_z = slicer.positions[slicer.triangles[id_last_contact_faces[j]][2]][2];
			double na = (v2.m_y - v1.m_y) * (v3.m_z - v1.m_z) - (v2.m_z - v1.m_z) * (v3.m_y - v1.m_y);
			double nb = (v2.m_z - v1.m_z) * (v3.m_x - v1.m_x) - (v2.m_x - v1.m_x) * (v3.m_z - v1.m_z);
			double nc = (v2.m_x - v1.m_x) * (v3.m_y - v1.m_y) - (v2.m_y - v1.m_y) * (v3.m_x - v1.m_x);
			Vector3 vn(na, nb, nc);
			vn.Normalized();
			last_contact_faces_normals.push_back(vn);
			// cout << "%%%%%%%%%%%" << last_contact_faces_normals[j].m_x << " " << last_contact_faces_normals[j].m_y << last_contact_faces_normals[j].m_z;
		}
		// cout << "%%%%%%%%%%%" << last_contact_faces_normals.size();
		cout << "cc" << endl;
		for (int i = 0; i < remove_triangles.size(); i++)
		{
			if (slicer.positions[remove_triangles[i][0]][2] != slicer.positions[remove_triangles[i][1]][2] || slicer.positions[remove_triangles[i][1]][2] != slicer.positions[remove_triangles[i][2]][2] || slicer.positions[remove_triangles[i][0]][2] != slicer.positions[remove_triangles[i][2]][2])
			{
				bool jud_contact = false;
				Vector3 v1, v2, v3;
				v1.m_x = slicer.positions[remove_triangles[i][0]][0];
				v1.m_y = slicer.positions[remove_triangles[i][0]][1];
				v1.m_z = slicer.positions[remove_triangles[i][0]][2];
				v2.m_x = slicer.positions[remove_triangles[i][1]][0];
				v2.m_y = slicer.positions[remove_triangles[i][1]][1];
				v2.m_z = slicer.positions[remove_triangles[i][1]][2];
				v3.m_x = slicer.positions[remove_triangles[i][2]][0];
				v3.m_y = slicer.positions[remove_triangles[i][2]][1];
				v3.m_z = slicer.positions[remove_triangles[i][2]][2];
				double ans = (v2.m_x - v1.m_x) * (v2.m_y - v3.m_y) - (v2.m_y - v1.m_y) * (v2.m_x - v3.m_x);
				if (ans > 0) // is clockwise
					swap(v2, v3);
				double na = (v2.m_y - v1.m_y) * (v3.m_z - v1.m_z) - (v2.m_z - v1.m_z) * (v3.m_y - v1.m_y);
				double nb = (v2.m_z - v1.m_z) * (v3.m_x - v1.m_x) - (v2.m_x - v1.m_x) * (v3.m_z - v1.m_z);
				double nc = (v2.m_x - v1.m_x) * (v3.m_y - v1.m_y) - (v2.m_y - v1.m_y) * (v3.m_x - v1.m_x);
				Vector3 vn(na, nb, nc);
				vn.Normalized();
				for (int j = 0; j < num_contact_faces; j++)
				{
					if (abs(vn.m_x - last_contact_faces_normals[j].m_x) < 0.001 && abs(vn.m_y - last_contact_faces_normals[j].m_y) < 0.001 && abs(vn.m_z - last_contact_faces_normals[j].m_z) < 0.001)
						jud_contact = true;
				}
				if (jud_contact == false)
					current_remove_triangles.push_back(remove_triangles[i]);
			}
		}
	}
	else
	{
		for (int i = 0; i < remove_triangles.size(); i++)
		{
			if (slicer.positions[remove_triangles[i][0]][2] != slicer.positions[remove_triangles[i][1]][2] || slicer.positions[remove_triangles[i][1]][2] != slicer.positions[remove_triangles[i][2]][2] || slicer.positions[remove_triangles[i][0]][2] != slicer.positions[remove_triangles[i][2]][2])
			{
				current_remove_triangles.push_back(remove_triangles[i]);
			}
		}
	}
	cout << "ccc" << endl;

	rotMatrix = Eigen::Quaterniond::FromTwoVectors(vector_after, vectorBefore).toRotationMatrix();
	for (int i = 0; i < slicer.positions.size(); i++)
	{
		Eigen::MatrixXd temp_V;
		temp_V.resize(3, 1);
		temp_V(0, 0) = slicer.positions[i][0];
		temp_V(1, 0) = slicer.positions[i][1];
		temp_V(2, 0) = slicer.positions[i][2];
		temp_V = rotMatrix.inverse() * temp_V;
		slicer.positions[i][0] = temp_V(0, 0);
		slicer.positions[i][1] = temp_V(1, 0);
		slicer.positions[i][2] = temp_V(2, 0);
	}
	current_slicer = slicer;

	// cout << current_remove_triangles.size() << endl;
	// if (height_of_beam_search != 2 || id_continue != 1)
	for (int i = 0; i < remove_triangles.size(); i++)
	{
		for (int j = 0; j < all_slicer.triangles.size();)
		{
			if (remove_triangles[i][0] == all_slicer.triangles[j][0] && remove_triangles[i][1] == all_slicer.triangles[j][1] && remove_triangles[i][2] == all_slicer.triangles[j][2])
			{
				all_slicer.triangles.erase(all_slicer.triangles.begin() + j);
				break;
			}
			j++;
		}
	}

	// add cutting plane triangles
	/*all_slicer.triangles.insert(all_slicer.triangles.begin(),cutting_plane_points)
	cutting_plane_points*/

	// erase some cutting plane
	for (int i = 0; i < all_cut_layers.size(); i++)
	{
		if (all_cut_layers_dependency_layer[i] == 0)
		{
			all_cut_layers.erase(all_cut_layers.begin() + i);
			all_cut_layers_dependency_layer.erase(all_cut_layers_dependency_layer.begin() + i);
			cutting_plane_edges.erase(cutting_plane_edges.begin() + i);
			cutting_plane_points.erase(cutting_plane_points.begin() + i);
			i--;
		}
	}
	cout << "cccc" << endl;
	// sort cutting_plane_points by adjacency relation
	vector<vector<int>> real_cutting_plane_triangles;
	real_cutting_plane_triangles.resize(all_cut_layers.size());
	for (int i = 0; i < all_cut_layers.size(); i++)
	{
		cv::Point2d current_triangle_point;
		cv::Point2d current_layer_point;
		current_layer_point.x = all_cut_layers[i][0].x;
		current_layer_point.y = all_cut_layers[i][0].y;
		double min_dis = MAX_D;
		int index_start_point_id;
		for (int j = 0; j < cutting_plane_edges[i].size(); j++)
		{
			current_triangle_point.x = all_slicer.positions[cutting_plane_edges[i][j].second][0];
			current_triangle_point.y = all_slicer.positions[cutting_plane_edges[i][j].second][1];
			if (distance_2d(current_layer_point, current_triangle_point) < min_dis)
			{
				index_start_point_id = j;
				min_dis = distance_2d(current_layer_point, current_triangle_point);
			}
		}
		int index_start_point = cutting_plane_edges[i][index_start_point_id].second;
		int index_of_last_edge = index_start_point_id;
		real_cutting_plane_triangles[i].push_back(cutting_plane_edges[i][index_start_point_id].first);
		// cout << "kkkk" << cutting_plane_edges[i].size() << endl;
		int jud_select_id = -1;
		int cont_segment = 0;
		while (index_start_point != cutting_plane_edges[i][index_start_point_id].first)
		{
			bool jud_select = false;
			cont_segment++;
			for (int j = 0; j < cutting_plane_edges[i].size(); j++)
			{
				if (j != index_of_last_edge && cutting_plane_edges[i][j].first == index_start_point)
				{
					real_cutting_plane_triangles[i].push_back(cutting_plane_edges[i][j].first);
					index_start_point = cutting_plane_edges[i][j].second;
					index_of_last_edge = j;
					jud_select = true;
					jud_select_id = 0;
					break;
				}
				else if (j != index_of_last_edge && cutting_plane_edges[i][j].second == index_start_point)
				{
					real_cutting_plane_triangles[i].push_back(cutting_plane_edges[i][j].second);
					index_start_point = cutting_plane_edges[i][j].first;
					index_of_last_edge = j;
					jud_select = true;
					jud_select_id = 1;
					break;
				}
			}
			if (cont_segment > 1000000)
			{
				jud_select = false;
				// real_cutting_plane_triangles[i].clear();
			}
			if (jud_select == false)
			{
				// cout << "no:" << i << "  ";
				jud_error[id_node] = true;
				/*for (int j = 0; j < cutting_plane_edges[i].size(); j++) {
					if (j != index_of_last_edge) {
						cv::Point2d current_point_2, current_point_3, current_point_4;
						current_point_3.x = all_slicer.positions[cutting_plane_edges[i][j].first][0];
						current_point_3.y = all_slicer.positions[cutting_plane_edges[i][j].first][1];
						current_point_4.x = all_slicer.positions[cutting_plane_edges[i][j].second][0];
						current_point_4.y = all_slicer.positions[cutting_plane_edges[i][j].second][1];
						if (jud_select_id == 0) {
							current_point_2.x = all_slicer.positions[cutting_plane_edges[i][index_of_last_edge].second][0];
							current_point_2.y = all_slicer.positions[cutting_plane_edges[i][index_of_last_edge].second][1];
						}
						else if (jud_select_id == 1) {
							current_point_2.x = all_slicer.positions[cutting_plane_edges[i][index_of_last_edge].first][0];
							current_point_2.y = all_slicer.positions[cutting_plane_edges[i][index_of_last_edge].first][1];
						}
						if (distance_2d(current_point_2, current_point_3) < 0.00001) {
							real_cutting_plane_triangles[i].push_back(cutting_plane_edges[i][j].first);
							index_start_point = cutting_plane_edges[i][j].second;
							index_of_last_edge = j;
							jud_select = true;
							jud_select_id = 0;
							break;
						}
						else if (distance_2d(current_point_2, current_point_4) < 0.00001) {
							real_cutting_plane_triangles[i].push_back(cutting_plane_edges[i][j].second);
							index_start_point = cutting_plane_edges[i][j].first;
							index_of_last_edge = j;
							jud_select = true;
							jud_select_id = 1;
							break;
						}
					}
				}*/
				break;
			}
		}
	}
	/*if (jud_error[id_node] == true)
		return;*/
	/*cout << "mm" << endl;*/
	/*if (height_of_beam_search == 6)
		return;*/
	for (int i = 0; i < real_cutting_plane_triangles.size(); i++)
	{
		if (real_cutting_plane_triangles[i].size() == 1)
		{
			real_cutting_plane_triangles.erase(real_cutting_plane_triangles.begin() + i);
			i--;
		}
	}

	///////////////////////////terminate//////////////////////////////
	if (Check_terminate(all_slicer))
	{
		remove_triangles = slicer.triangles;
		jud_outer_beam_search_terminate = true;
	}

	//////////////////////////////////////////////////////////////////

	vector<int> id_contact_faces;
	// if (height_of_beam_search != 2 || id_continue != 1)
	if (using_solid_model == true)
	{
		Anticlockwise(real_cutting_plane_triangles, all_slicer);
		using Coord = double;
		using NN = uint32_t;
		using Point = std::array<Coord, 2>;
		for (int i = 0; i < real_cutting_plane_triangles.size(); i++)
		{
			if (flag_cut_layers_is_hole[i] != -1)
				continue;
			std::vector<std::vector<Point>> polygon;
			polygon.clear();
			std::vector<Point> temp_vec;
			polygon.push_back(temp_vec);
			map<int, int> map_index_faces;
			map_index_faces.clear();
			for (int j = 0; j < real_cutting_plane_triangles[i].size(); j++)
			{
				polygon[0].push_back({all_slicer.positions[real_cutting_plane_triangles[i][j]][0], all_slicer.positions[real_cutting_plane_triangles[i][j]][1]});
				map_index_faces.insert({j, real_cutting_plane_triangles[i][j]});
			}
			for (int m = 0; m < real_cutting_plane_triangles.size(); m++)
			{
				if (flag_cut_layers_is_hole[m] != -1 && follow_index[flag_cut_layers_is_hole[m]] == i)
				{
					for (int j = 0; j < real_cutting_plane_triangles[m].size(); j++)
					{
						polygon.push_back(temp_vec);
						polygon[1].push_back({all_slicer.positions[real_cutting_plane_triangles[m][j]][0], all_slicer.positions[real_cutting_plane_triangles[m][j]][1]});
						map_index_faces.insert({polygon[0].size() + j, real_cutting_plane_triangles[m][j]});
					}
					// 多个contour?
					break;
				}
			}

			std::vector<NN> indices = mapbox::earcut<NN>(polygon);

			for (int j = 0; j < indices.size();)
			{
				TRiangle the_new_cutting_plane_triangle;
				the_new_cutting_plane_triangle[0] = map_index_faces[indices[j]];
				j++;
				the_new_cutting_plane_triangle[1] = map_index_faces[indices[j]];
				j++;
				the_new_cutting_plane_triangle[2] = map_index_faces[indices[j]];
				j++;
				all_slicer.triangles.insert(all_slicer.triangles.end(), the_new_cutting_plane_triangle); // 添加接触面
				id_contact_faces.push_back(all_slicer.triangles.size() - 1);
			}
		}
		//}
	}

	/*if (height_of_beam_search == 1) {
		for (int i = 0; i < real_cutting_plane_triangles.size(); i++) {
			VEctor new_point;
			double temp_x = 0, temp_y = 0;
			for (int j = 0; j < real_cutting_plane_triangles[i].size(); j++) {
				temp_x += all_slicer.positions[real_cutting_plane_triangles[i][j]][0];
				temp_y += all_slicer.positions[real_cutting_plane_triangles[i][j]][1];
			}
			new_point[0] = temp_x / real_cutting_plane_triangles[i].size();
			new_point[1] = temp_y / real_cutting_plane_triangles[i].size();
			new_point[2] = all_slicer.positions[real_cutting_plane_triangles[i][0]][2];
			all_slicer.positions.insert(all_slicer.positions.end(), new_point);
			for (int j = 0; j < real_cutting_plane_triangles[i].size(); j++) {
				TRiangle the_new_cutting_plane_triangle;
				the_new_cutting_plane_triangle[0] = all_slicer.positions.size() - 1;
				the_new_cutting_plane_triangle[1] = real_cutting_plane_triangles[i][j];
				the_new_cutting_plane_triangle[2] = real_cutting_plane_triangles[i][(j + 1) % real_cutting_plane_triangles[i].size()];
				all_slicer.triangles.insert(all_slicer.triangles.end(), the_new_cutting_plane_triangle);
			}
		}
	}*/

	rotMatrix = Eigen::Quaterniond::FromTwoVectors(vector_after, vectorBefore).toRotationMatrix();
	//////rotate back //////
	// if (height_of_beam_search != 2)
	for (int i = 0; i < all_slicer.positions.size(); i++)
	{
		Eigen::MatrixXd temp_V;
		temp_V.resize(3, 1);
		temp_V(0, 0) = all_slicer.positions[i][0];

		temp_V(1, 0) = all_slicer.positions[i][1];
		temp_V(2, 0) = all_slicer.positions[i][2];
		temp_V = rotMatrix.inverse() * temp_V;
		all_slicer.positions[i][0] = temp_V(0, 0);
		all_slicer.positions[i][1] = temp_V(1, 0);
		all_slicer.positions[i][2] = temp_V(2, 0);
	}
	/////////////////////////

	Slicer_2 all_slicer_2 = all_slicer;
	all_slicer_2.triangles = remove_triangles;
	if (judge_continue_additive == false)
		all_slicer_2.save(file_name + "-" + to_string(height_of_beam_search) + "_" + to_string(cont_number_of_queue) + "_current" + ".obj");
	else
		all_slicer_2.save(file_name + "-" + to_string(height_of_beam_search) + "_" + to_string(cont_number_of_queue) + "_" + to_string(id_continue) + "_current" + "_subblock.obj");

	Eigen::MatrixXd temp_V(all_slicer.positions.size(), 3);
	Eigen::MatrixXi temp_F(all_slicer.triangles.size(), 3);
	for (int i = 0; i < all_slicer.positions.size(); i++)
	{
		temp_V.row(i).x() = all_slicer.positions[i][0];
		temp_V.row(i).y() = all_slicer.positions[i][1];
		temp_V.row(i).z() = all_slicer.positions[i][2];
	}
	for (int i = 0; i < all_slicer.triangles.size(); i++)
	{
		temp_F.row(i).x() = all_slicer.triangles[i][0];
		temp_F.row(i).y() = all_slicer.triangles[i][1];
		temp_F.row(i).z() = all_slicer.triangles[i][2];
	}
	if (judge_continue_additive == false)
	{
		Katana::Instance().stl.saveStlFromObj(file_name + "-" + to_string(height_of_beam_search) + "_" + to_string(cont_number_of_queue) + ".stl", temp_V, temp_F);
		vector<std::array<double, 3>> V3;
		vector<std::array<int, 3>> F3;
		vector<std::array<double, 3>> N3;
		ifstream ifs(file_name + "-" + to_string(height_of_beam_search) + "_" + to_string(cont_number_of_queue) + ".stl");
		igl::read_stl_ascii(ifs, V3, F3, N3);
		Eigen::MatrixXd V4(V3.size(), 3);
		Eigen::MatrixXi F4(F3.size(), 3);
		Eigen::MatrixXd N4(N3.size(), 3);
		for (int i = 0; i < V3.size(); i++)
			for (int j = 0; j < 3; j++)
				V4(i, j) = V3[i][j];
		for (int i = 0; i < F3.size(); i++)
			for (int j = 0; j < 3; j++)
				F4(i, j) = F3[i][j];
		igl::writeSTL(file_name + "-" + to_string(height_of_beam_search) + "_" + to_string(cont_number_of_queue) + "_B.stl", V4, F4, igl::FileEncoding::Binary);

		Geometry tessel;
		tessel.visit(ImportSTL(file_name + "-" + to_string(height_of_beam_search) + "_" + to_string(cont_number_of_queue) + "_B.stl"));
		tessel.visit(ExportOBJ(file_name + "-" + to_string(height_of_beam_search) + "_" + to_string(cont_number_of_queue) + ".obj"));
	}
	else
	{
		Katana::Instance().stl.saveStlFromObj(file_name + "-" + to_string(height_of_beam_search) + "_" + to_string(cont_number_of_queue) + "_" + to_string(id_continue) + "_subblock.stl", temp_V, temp_F);
		vector<std::array<double, 3>> V3;
		vector<std::array<int, 3>> F3;
		vector<std::array<double, 3>> N3;
		ifstream ifs(file_name + "-" + to_string(height_of_beam_search) + "_" + to_string(cont_number_of_queue) + "_" + to_string(id_continue) + "_subblock.stl");
		igl::read_stl_ascii(ifs, V3, F3, N3);
		Eigen::MatrixXd V4(V3.size(), 3);
		Eigen::MatrixXi F4(F3.size(), 3);
		Eigen::MatrixXd N4(N3.size(), 3);
		for (int i = 0; i < V3.size(); i++)
			for (int j = 0; j < 3; j++)
				V4(i, j) = V3[i][j];
		for (int i = 0; i < F3.size(); i++)
			for (int j = 0; j < 3; j++)
				F4(i, j) = F3[i][j];
		igl::writeSTL(file_name + "-" + to_string(height_of_beam_search) + "_" + to_string(cont_number_of_queue) + "_" + to_string(id_continue) + "_subblock_B.stl", V4, F4, igl::FileEncoding::Binary);

		Geometry tessel;
		tessel.visit(ImportSTL(file_name + "-" + to_string(height_of_beam_search) + "_" + to_string(cont_number_of_queue) + "_" + to_string(id_continue) + "_subblock_B.stl"));
		tessel.visit(ExportOBJ(file_name + "-" + to_string(height_of_beam_search) + "_" + to_string(cont_number_of_queue) + "_" + to_string(id_continue) + "_subblock.obj"));
	}

	string str_contact_faces;
	if (judge_continue_additive == false)
		str_contact_faces = file_name + "-" + to_string(height_of_beam_search) + "_" + to_string(cont_number_of_queue) + "_contact_faces.txt";
	else
		str_contact_faces = file_name + "-" + to_string(height_of_beam_search) + "_" + to_string(cont_number_of_queue) + "_" + to_string(id_continue) + "_subblock_contact_faces.txt";

	Visual vis_2;
	cv::Point3d input_ori(vector_after.x(), vector_after.y(), vector_after.z());
	if (judge_continue_additive == false)
		vis_2.generateArrows(input_ori, file_name + "-" + to_string(height_of_beam_search) + "_" + to_string(cont_number_of_queue));
	else
		vis_2.generateArrows(input_ori, file_name + "-" + to_string(height_of_beam_search) + "_" + to_string(cont_number_of_queue) + "_" + to_string(id_continue) + "_subblock");

	ofstream ofile_contact_faces(str_contact_faces);
	ofile_contact_faces << id_contact_faces.size() << endl;
	for (int i = 0; i < id_contact_faces.size(); i++)
		ofile_contact_faces << id_contact_faces[i] << endl;
	cout << "d" << endl;
	/////////////////////////
	return;
}

vector<vector<int>> HybridManufacturing::FindAllCutLayers(Layer_Graph layer_graph, vector<vector<int>> final_pathes, vector<vector<int>> &all_cut_layers_dependency_layer, bool &jud_admit)
{
	vector<vector<int>> all_id_of_cut_layers;
	all_id_of_cut_layers.resize(final_pathes.size());
	all_cut_layers_dependency_layer.resize(final_pathes.size());
	for (int i = 0; i < final_pathes.size(); i++)
	{
		for (int j = 0; j < final_pathes[i].size(); j++)
		{
			bool jud_cut_layer = true;
			bool jud_cut_layer_fork = false;
			for (int k = 0; k < layer_graph.G_2[final_pathes[i][j]].size(); k++)
			{
				bool jud_cut_layer_2 = true;
				for (int m = 0; m < final_pathes[i].size(); m++)
				{
					if (j != m)
					{
						if (layer_graph.edges[layer_graph.G_2[final_pathes[i][j]][k]].GetFrom() == final_pathes[i][m])
						{
							jud_cut_layer_2 = false;
							break;
						}
					}
				}
				if (jud_cut_layer_2 == false)
				{
					jud_cut_layer = false;
				}
				if (jud_cut_layer_2 == true)
				{
					jud_cut_layer_fork = true;
				}
			}
			if (jud_cut_layer == false && jud_cut_layer_fork == true)
			{ // 对于中间存在分叉的情况，先剔除这些结果
				jud_admit = false;
				return all_id_of_cut_layers;
			}
			if (jud_cut_layer == true || layer_graph.G_2[final_pathes[i][j]].size() == 0)
			{
				if (layer_graph.G_2[final_pathes[i][j]].size() == 0 && layer_graph.G[final_pathes[i][j]].size() == 0) // 删除孤立layer
					continue;
				all_id_of_cut_layers[i].push_back(j);
				all_cut_layers_dependency_layer[i].push_back(layer_graph.G_2[final_pathes[i][j]].size());
			}
		}
	}
	return all_id_of_cut_layers;
}

void HybridManufacturing::sort_candidate_nodes(vector<int> &candidate_nodes, vector<vector<int>> Tree_nodes_for_S)
{
	/*for (int i = 0; i < candidate_nodes.size(); i++) {
		for (int j = i + 1; j < candidate_nodes.size(); j++) {
			if (Tree_nodes_for_S[candidate_nodes[i]].size() < Tree_nodes_for_S[candidate_nodes[j]].size()) {
				swap(candidate_nodes[i], candidate_nodes[j]);
			}
		}
	}*/
}

void HybridManufacturing::outputToObj(vector<Point3ff> outputVertices, vector<face> outputFace, Point3ff outputColor, string filename)
{
	// 检查是否需要输出颜色
	ofstream fout;
	fout.open(filename);
	if (outputColor == Point3ff(-1, -1, -1))
	{
		for (int i = 0; i < outputVertices.size(); i++)
		{
			fout << "v " << outputVertices[i].x << " " << outputVertices[i].y << " " << outputVertices[i].z << endl;
		}
		for (int i = 0; i < outputFace.size(); i++)
		{
			fout << "f " << outputFace[i].v1 + 1 << " " << outputFace[i].v2 + 1 << " " << outputFace[i].v3 + 1 << endl;
		}
	}
	else
	{
		for (int i = 0; i < outputVertices.size(); i++)
		{
			fout << "v " << outputVertices[i].x << " " << outputVertices[i].y << " " << outputVertices[i].z << " ";
			fout << outputColor.x << " " << outputColor.y << " " << outputColor.z << endl;
		}
		for (int i = 0; i < outputFace.size(); i++)
		{
			fout << "f " << outputFace[i].v1 + 1 << " " << outputFace[i].v2 + 1 << " " << outputFace[i].v3 + 1 << endl;
		}
	}
}

// 判断两个向量几乎相同
bool isEqual(Point3ff v1, Point3ff v2)
{
	if (v1.x - v2.x > -1e-2 && v1.x - v2.x < 1e-2)
		if (v1.y - v2.y > -1e-2 && v1.y - v2.y < 1e-2)
			if (v1.z - v2.z > -1e-2 && v1.z - v2.z < 1e-2)
				return true;
	return false;
}

// 判断两个向量方向相同
// 因为经常会有法向跟打印方向差了0.0几导致不相等,所以写个这个
bool isSameDirection(Point3ff v1, Point3ff v2)
{
	if (Dot(v1, v2) > 0)
		return true;
	else
		return false;
}

// 在-direction方向下看去,三点是否为逆时针
int isAnticlockwise(Point3ff direction, Point3ff p1, Point3ff p2, Point3ff p3)
{
	// p2p1 叉乘 p2p3 与direction方向相同 则为逆时针
	Point3ff p21 = p1 - p2;
	Point3ff p23 = p3 - p2;
	// 判断是否为一条直线
	p21 = p21 / length_P(p21, Point3ff(0, 0, 0));
	p23 = p23 / length_P(p23, Point3ff(0, 0, 0));
	Point3ff ans = Cross(p21, p23);

	if (isSameDirection(ans, direction * (-1))) // is anticlockwise
		return 1;
	else if (isEqual(ans, Point3ff(0, 0, 0)) || isEqual(p21 + p23, Point3ff(0, 0, 0)))
		return 2;
	// ans=sin(p21,p23) sin1度为0.017
	else if (length_P(ans, Point3ff(0, 0, 0)) > -0.02 && length_P(ans, Point3ff(0, 0, 0)) < 0.02)
		return 2;
	else // is clockwise
		return 3;
}

// 检测三点中间是否还有其他顶点
bool checkHaveNoOtherPoint(Point3ff direction, Point3ff p1, Point3ff p2, Point3ff p3, std::list<Point3ff> pointList)
{
	std::list<Point3ff>::iterator iter = pointList.begin();
	for (; iter != pointList.end(); iter++)
	{
		Point3ff p = *iter;
		if (p == p1 || p == p2 || p == p3)
			continue;
		// if point p is in the triangle, (p,p1,p2),(p,p2,p3),(p,p3,p1) should be anticlockwise
		if (isAnticlockwise(direction, p, p1, p2) == 1 && isAnticlockwise(direction, p, p2, p3) == 1 && isAnticlockwise(direction, p, p3, p1) == 1)
			return false;
	}
	return true;
}

// 计算三角形有向面积
double calculateTriArea(Point3ff direction, Point3ff p1, Point3ff p2, Point3ff p3)
{
	Point3ff N = Cross(p2 - p1, p3 - p1);
	double S = 0.5 * length_P(N, Point3ff(0, 0, 0));
	if (isSameDirection(N, direction))
		return S;
	else
		return -S;
}

// 判断这一层顶点在-direction方向下看去是否为逆时针
bool pointsAreAnticlockwise(Point3ff direction, vector<Point3ff> inputVertices)
{
	// 输入检查
	if (inputVertices.size() < 3)
		return false;
	// 计算有向面积
	double S = 0;
	Point3ff p1 = inputVertices[0];
	for (int i = 1; i <= inputVertices.size() - 3; i++)
	{
		Point3ff p2 = inputVertices[i];
		Point3ff p3 = inputVertices[i + 1];
		S = S + calculateTriArea(direction, p1, p2, p3);
	}
	if (S > 0)
		return true;
	else
		return false;
}

// 任意(凹)多边形转换为多个三角形
// 顶面在-direction方向下看去是逆时针
// 底面在direction方向下看去是逆时针
// inputVertices此时为一层多边形轮廓
void polygonToTriangle(Point3ff direction, vector<Point3ff> inputVertices, bool isAnti, double h, vector<Point3ff> &outputVertices, vector<face> &outputFace)
{
	int numOfPoint = inputVertices.size() - 1;
	// 输入检查
	if (numOfPoint < 3)
		return;
	// 放入列表,应在direction方向下加h层厚(顶面+h,底面h=0)
	list<Point3ff> pointList;
	if (!isAnti)
	{
		for (int i = numOfPoint - 1; i >= 0; i--)
			pointList.push_back(inputVertices[i] + direction * h);
	}
	else
	{
		for (int i = 0; i < numOfPoint; i++)
			pointList.push_back(inputVertices[i] + direction * h);
	}
	// 一些2021.12.2写的函数
	std::list<Point3ff>::iterator iter = pointList.begin();
	Point3ff p1, p2, p3;
	p1 = *(iter);
	p2 = *(next(iter));
	p3 = *(next(iter, 2));
	while (pointList.size() >= 3)
	{
		// 最后三个顶点特殊处理
		if (pointList.size() == 3)
		{
			if (isAnticlockwise(direction, p1, p2, p3) == 1)
			{
				// use face (p1,p2,p3)
				int indexOfP = outputVertices.size();
				outputVertices.push_back(p1);
				outputVertices.push_back(p2);
				outputVertices.push_back(p3);
				face theFace(indexOfP, indexOfP + 1, indexOfP + 2);
				outputFace.push_back(theFace);
			}
			else if (isAnticlockwise(direction, p1, p2, p3) == 3)
			{
				// use face (p1,p3,p2)
				int indexOfP = outputVertices.size();
				outputVertices.push_back(p1);
				outputVertices.push_back(p3);
				outputVertices.push_back(p2);
				face theFace(indexOfP, indexOfP + 1, indexOfP + 2);
				outputFace.push_back(theFace);
			}
			else // 不大应该出现
			{
				cout << "最后三个顶点是一条线.请特殊处理." << endl;
			}
			break;
		}
		// if(p1,p2,p3) is anticlockwise and no other point in this polygon
		if (isAnticlockwise(direction, p1, p2, p3) == 1 && checkHaveNoOtherPoint(direction, p1, p2, p3, pointList))
		{
			// use face (p1,p2,p3)
			int indexOfP = outputVertices.size();
			outputVertices.push_back(p1);
			outputVertices.push_back(p2);
			outputVertices.push_back(p3);
			face theFace(indexOfP, indexOfP + 1, indexOfP + 2);
			outputFace.push_back(theFace);
			// delete p2
			if (next(iter) == pointList.end())
				pointList.erase(pointList.begin());
			else
				pointList.erase(next(iter));
			if (pointList.size() >= 3)
			{
				if (next(iter) == pointList.end())
				{
					p2 = *(pointList.begin());
					p3 = *(++pointList.begin());
				}
				else if (next(iter, 2) == pointList.end())
				{
					p2 = *(next(iter));
					p3 = *(pointList.begin());
				}
				else
				{
					p2 = *(next(iter));
					p3 = *(next(iter, 2));
				}
			}
		}
		else if (isAnticlockwise(direction, p1, p2, p3) == 2)
		{
			// p1->p2->p3 is a line, so delete p2
			if (next(iter) == pointList.end())
				pointList.erase(pointList.begin());
			else
				pointList.erase(next(iter));
			if (pointList.size() >= 3)
			{
				if (next(iter) == pointList.end())
				{
					p2 = *(pointList.begin());
					p3 = *(++pointList.begin());
				}
				else if (next(iter, 2) == pointList.end())
				{
					p2 = *(next(iter));
					p3 = *(pointList.begin());
				}
				else
				{
					p2 = *(next(iter));
					p3 = *(next(iter, 2));
				}
			}
		}
		else
		{
			iter++;
			if (iter == pointList.end())
			{
				iter = pointList.begin();
				p1 = *(iter);
				p2 = *(next(iter));
				p3 = *(next(iter, 2));
			}
			else if (next(iter) == pointList.end())
			{
				p1 = *(iter);
				p2 = *(pointList.begin());
				p3 = *(++pointList.begin());
			}
			else if (next(iter, 2) == pointList.end())
			{
				p1 = *(iter);
				p2 = *(next(iter));
				p3 = *(pointList.begin());
			}
			else
			{
				p1 = *(iter);
				p2 = *(next(iter));
				p3 = *(next(iter, 2));
			}
		}
	}
}

// 生成侧壁
void generateWall(Point3ff direction, vector<Point3ff> inputVertices, bool isAnti, double h, vector<Point3ff> &outputVertices, vector<face> &outputFace)
{
	// 对每一对顶点x和y,生成打印方向对应的xx和yy,然后连接成两个三角形面片
	for (int xIndex = 0; xIndex < inputVertices.size() - 1; xIndex++)
	{
		// 构造顶点
		Point3ff x = inputVertices[xIndex];
		Point3ff y = inputVertices[xIndex + 1];
		Point3ff xx = x + direction * h;
		Point3ff yy = y + direction * h;
		int indexOfX = outputVertices.size();
		outputVertices.push_back(x);  // indexOfX
		outputVertices.push_back(xx); // indexOfX+1
		outputVertices.push_back(y);  // indexOfX+2
		outputVertices.push_back(yy); // indexOfX+3

		if (isAnti)
		{
			// 构造面片[x,y,xx]和[y,yy,xx]
			face f1(indexOfX, indexOfX + 2, indexOfX + 1);
			face f2(indexOfX + 2, indexOfX + 3, indexOfX + 1);
			outputFace.push_back(f1);
			outputFace.push_back(f2);
		}
		else
		{
			// 构造面片[x,xx,y]和[y,xx,yy]
			face f1(indexOfX, indexOfX + 1, indexOfX + 2);
			face f2(indexOfX + 2, indexOfX + 1, indexOfX + 3);
			outputFace.push_back(f1);
			outputFace.push_back(f2);
		}
	}
}

// 生成环
// 顶面在-direction方向下看去是逆时针
// 底面在direction方向下看去是逆时针
void generateLoop(Point3ff direction, vector<Point3ff> contour1, vector<Point3ff> contour2, double h, vector<Point3ff> &outputVertices, vector<face> &outputFace)
{
	// 这个函数里会默认数组首尾顶点是不同的,所以需要删掉最后一个顶点
	contour1.pop_back();
	contour2.pop_back();
	// 确认起始点,按最近策略连接顶点
	double dis = length(contour1[0], contour2[0]);
	int indexOfVertex = 0; // 挑一个离contour1[0]最近的顶点
	for (int i = 1; i < contour2.size(); i++)
	{
		if (dis > length(contour1[0], contour2[i]))
		{
			dis = length(contour1[0], contour2[i]);
			indexOfVertex = i;
		}
	}
	// 重新排序contour2,此时两个起始点最近
	vector<Point3ff> temp;
	for (int i = 0; i < contour2.size(); i++)
		temp.push_back(contour2[(i + indexOfVertex) % contour2.size()]);
	contour2 = temp;

	// 循环:每次拿出一个contour1和一个contour2的顶点,并在contour1[i+1]和contour2[j+1]这两个候选顶点中选择一个,连接为三角形
	// 选择方法:数目差/组成的三角形内部不包括候选顶点和环中心/三点不在一条直线上
	// 退出条件:有个三角形的顶点是初始两个顶点之一,即连接了一圈回来了
	int index1 = 0, index2 = 0;
	while (index1 != contour1.size() && index2 != contour2.size())
	{
		Point3ff p1 = contour1[index1];
		Point3ff p2 = contour2[index2];
		Point3ff p3 = contour1[(index1 + 1) % contour1.size()];
		bool chooseContour1 = true; // p3选contour1上的
		// 首先比较剩余顶点的数目,优先选剩的多的
		int numOfV1 = contour1.size() - index1, numOfV2 = contour2.size() - index2;
		if (numOfV1 < numOfV2)
		{
			p3 = contour2[(index2 + 1) % contour2.size()];
			chooseContour1 = false;
		}
		// 接下来判断该三角形内部不包含候选顶点和环中心
		bool check = true;

		if (chooseContour1) // 外-内-外应该是顺时针,且不包含contour2[index2+1],且该点和p1p3,p2p3不在一条直线上
		{
			if (isAnticlockwise(direction, p1, p2, p3) != 3)
				check = false;
			Point3ff p = contour2[(index2 + 1) % contour2.size()];
			if (isAnticlockwise(direction, p, p2, p1) == 1 && isAnticlockwise(direction, p, p3, p2) == 1 && isAnticlockwise(direction, p, p1, p3) == 1)
				check = false;
			if (isAnticlockwise(direction, p1, p, p3) == 2 || isAnticlockwise(direction, p2, p, p3) == 2)
				check = false;
		}
		else // 外-内-内应该是顺时针,且不包含contour1[index1+1],且该点和p1p3,p2p3不在一条直线上
		{
			if (isAnticlockwise(direction, p1, p2, p3) != 3)
				check = false;
			Point3ff p = contour1[(index1 + 1) % contour1.size()];
			if (isAnticlockwise(direction, p, p2, p1) == 1 && isAnticlockwise(direction, p, p3, p2) == 1 && isAnticlockwise(direction, p, p1, p3) == 1)
				check = false;
			if (isAnticlockwise(direction, p1, p, p3) == 2 || isAnticlockwise(direction, p2, p, p3) == 2)
				check = false;
		}

		// 总能选出一个来吧,就不加两个都不能选的特殊处理了
		if (check == false && chooseContour1 == true)
		{
			chooseContour1 = false;
			p3 = contour2[(index2 + 1) % contour2.size()];
		}
		else if (check == false && chooseContour1 == false)
		{
			chooseContour1 = true;
			p3 = contour1[(index1 + 1) % contour1.size()];
		}

		// 顶面:将该面片[p1,p3,p2]+direction*h (逆时针)添加到数组中
		int indexOfP = outputVertices.size();
		outputVertices.push_back(p1 + direction * h);
		outputVertices.push_back(p3 + direction * h);
		outputVertices.push_back(p2 + direction * h);
		face theFace(indexOfP, indexOfP + 1, indexOfP + 2);
		outputFace.push_back(theFace);

		// 底面:将该面片[p1,p2,p3] (顺时针)添加到数组中
		indexOfP = outputVertices.size();
		outputVertices.push_back(p1);
		outputVertices.push_back(p2);
		outputVertices.push_back(p3);
		theFace = face(indexOfP, indexOfP + 1, indexOfP + 2);
		outputFace.push_back(theFace);

		// 这里不对数组大小取余
		if (chooseContour1)
			index1++;
		else
			index2++;
	}
	// 最后的特殊处理,此时有一侧contour连接完毕,将另一侧连接过来
	if (index1 == contour1.size())
	{
		while (index2 != contour2.size())
		{
			Point3ff p1 = contour1[0];
			Point3ff p2 = contour2[index2];
			Point3ff p3 = contour2[(index2 + 1) % contour2.size()];

			// 应该不需要判断能否加入吧

			// 顶面:将该面片[p1,p3,p2]+direction*h (逆时针)添加到数组中
			int indexOfP = outputVertices.size();
			outputVertices.push_back(p1 + direction * h);
			outputVertices.push_back(p3 + direction * h);
			outputVertices.push_back(p2 + direction * h);
			face theFace(indexOfP, indexOfP + 1, indexOfP + 2);
			outputFace.push_back(theFace);

			// 底面:将该面片[p1,p2,p3] (顺时针)添加到数组中
			indexOfP = outputVertices.size();
			outputVertices.push_back(p1);
			outputVertices.push_back(p2);
			outputVertices.push_back(p3);
			theFace = face(indexOfP, indexOfP + 1, indexOfP + 2);
			outputFace.push_back(theFace);

			index2++;
		}
	}
	else if (index2 == contour2.size())
	{
		while (index1 != contour1.size())
		{
			Point3ff p1 = contour1[index1];
			Point3ff p2 = contour2[0];
			Point3ff p3 = contour1[(index1 + 1) % contour1.size()];

			// 应该不需要判断能否加入吧

			// 顶面:将该面片[p1,p3,p2]+direction*h (逆时针)添加到数组中
			int indexOfP = outputVertices.size();
			outputVertices.push_back(p1 + direction * h);
			outputVertices.push_back(p3 + direction * h);
			outputVertices.push_back(p2 + direction * h);
			face theFace(indexOfP, indexOfP + 1, indexOfP + 2);
			outputFace.push_back(theFace);

			// 底面:将该面片[p1,p2,p3] (顺时针)添加到数组中
			indexOfP = outputVertices.size();
			outputVertices.push_back(p1);
			outputVertices.push_back(p2);
			outputVertices.push_back(p3);
			theFace = face(indexOfP, indexOfP + 1, indexOfP + 2);
			outputFace.push_back(theFace);

			index1++;
		}
	}
}

// 新-生成面片 支持环结构
void newGenerateMesh(Point3ff direction, bool isLoop, vector<Point3ff> contour1, vector<Point3ff> contour2, double h, vector<Point3ff> &outputVertices, vector<face> &outputFace)
{
	// isLoop为false,和原来一致,生成侧壁和顶面底面的三角面片
	// 此时contour2不使用
	if (!isLoop)
	{
		// 判断顶点方向
		bool isAnti = pointsAreAnticlockwise(direction, contour1);
		// 生成侧壁
		generateWall(direction, contour1, isAnti, h, outputVertices, outputFace);
		// 顶面传入h
		polygonToTriangle(direction, contour1, isAnti, h, outputVertices, outputFace);
		// 底面传入h=0
		polygonToTriangle(direction * (-1), contour1, !isAnti, 0, outputVertices, outputFace);
	}
	else // isLoop为true,为环结构,contour1为外层
	{
		// 有可能两层环方向不同,将其都变为在-direction方向下看去是逆时针
		if (!pointsAreAnticlockwise(direction, contour1))
		{
			vector<Point3ff> temp;
			for (int i = contour1.size() - 1; i >= 0; i--)
				temp.push_back(contour1[i]);
			contour1 = temp;
		}
		if (!pointsAreAnticlockwise(direction, contour2))
		{
			vector<Point3ff> temp;
			for (int i = contour2.size() - 1; i >= 0; i--)
				temp.push_back(contour2[i]);
			contour2 = temp;
		}
		bool isAnti = true;
		// 生成侧壁
		generateWall(direction, contour1, isAnti, h, outputVertices, outputFace);
		generateWall(direction, contour2, !isAnti, h, outputVertices, outputFace); // 注意法向是反的,所以输入了!isAnti
		// 生成顶面和地面
		generateLoop(direction, contour1, contour2, h, outputVertices, outputFace);
	}
}

// 生成侧壁和顶面底面的三角面片
void generateMesh(Point3ff direction, vector<vector<Point3ff>> inputVertices, double h, vector<Point3ff> &outputVertices, vector<face> &outputFace)
{
	// 对于layerIndex层,生成侧壁,底面和顶面
	for (int layerIndex = 0; layerIndex < inputVertices.size(); layerIndex++)
	{
		// 判断顶点方向
		bool isAnti = pointsAreAnticlockwise(direction, inputVertices[layerIndex]);
		// 生成侧壁
		generateWall(direction, inputVertices[layerIndex], isAnti, h, outputVertices, outputFace);
		// 顶面传入h
		polygonToTriangle(direction, inputVertices[layerIndex], isAnti, h, outputVertices, outputFace);
		// 底面传入h=0
		polygonToTriangle(direction * (-1), inputVertices[layerIndex], !isAnti, 0, outputVertices, outputFace);
		printf("正在处理%d/%d.\n", layerIndex + 1, (int)inputVertices.size());
	}
}

void HybridManufacturing::visualize_layers(string file_name, Point3ff vn, bool judge_continue_additive, int id_continue)
{
	cout << "aa" << endl;
	string input_file_name;
	if (judge_continue_additive == false)
		input_file_name = file_name + ".txt";
	else
		input_file_name = file_name + "_" + to_string(id_continue) + "_subblock.txt";
	cout << input_file_name << endl;
	ifstream ifile(input_file_name);
	int num_contours;
	ifile >> num_contours;
	vector<vector<Point3ff>> inputVertices;
	for (int i = 0; i < num_contours; i++)
	{
		vector<Point3ff> currentContour;
		int num_points;
		ifile >> num_points;
		for (int j = 0; j < num_points; j++)
		{
			// cout << "a";
			Point3ff currentP;
			ifile >> currentP.x >> currentP.y >> currentP.z;
			currentContour.push_back(currentP);
		}
		inputVertices.push_back(currentContour);
	}
	cout << "bb" << num_contours << endl;
	// 输入参数
	double R = 0.3; // 切向圆的半径
	int nr = 20;	// 用几边形来拟合横截面圆
	int type = 1;	// 划分成网格时的类型 0为四边形网格,暂不用
	//        b---d           b---d
	// mode 1 | / |    mode 2 | \ |  3为1和2交替 4为2和1交替
	//        a---c           a---c
	// 默认type=1

	double alpha = 1, beta = 1;		 // 椭圆,alpha*副法向B, beta*正法向N
	Point3ff outputColor(1, 1, 255); // 模型颜色,无颜色为(-1,-1,-1)
	string objPath;
	if (judge_continue_additive == false)
		objPath = file_name + "_vis_small.obj"; // 输出obj的路径
	else
		objPath = file_name + "_" + to_string(id_continue) + "_subblock_vis_small.obj"; // 输出obj的路径

	// 调用函数
	vector<Point3ff> outputVertices;
	vector<face> outputFace;
	// Point3ff v1(-10.886800, 8.700820, 64.546097), v2(-13.794400, 9.800380, 62.789299), v3(-22.245300, 2.629520, 68.324501);
	// double na = (v2.y - v1.y) * (v3.z - v1.z) - (v2.z - v1.z) * (v3.y - v1.y);
	// double nb = (v2.z - v1.z) * (v3.x - v1.x) - (v2.x - v1.x) * (v3.z - v1.z);
	// double nc = (v2.x - v1.x) * (v3.y - v1.y) - (v2.y - v1.y) * (v3.x - v1.x);
	// Point3ff vn(na, nb, nc);
	// double t = sqrt(vn.x * vn.x + vn.y * vn.y + vn.z * vn.z);
	// vn.x /= t;
	// vn.y /= t;
	// vn.z /= t;

	// Point3ff direction(0, 0, 1);	//单位向量
	generateClosedTube(vn, inputVertices, R, nr, type, alpha, beta, outputVertices, outputFace);
	cout << "cc" << endl;
	// 输出为obj
	outputToObj(outputVertices, outputFace, outputColor, objPath);
	return;
}

void HybridManufacturing::visualize_layers_stair_case(string file_name, string vis_file_contain, Point3ff vn, bool judge_continue_additive, int id_continue)
{
	string input_file_name, input_file_name_contain;
	if (judge_continue_additive == false)
		input_file_name = file_name + ".txt";
	else
		input_file_name = file_name + "_" + to_string(id_continue) + "_subblock.txt";

	// 输入参数
	double h = 2.0; // 每层的厚度

	double alpha = 1, beta = 1;		 // 椭圆,alpha*副法向B, beta*正法向N
	Point3ff outputColor(1, 1, 255); // 模型颜色,无颜色为(-1,-1,-1)
	string objPath;
	if (judge_continue_additive == false)
		objPath = file_name + "_vis_stair_case.obj"; // 输出obj的路径
	else
		objPath = file_name + "_" + to_string(id_continue) + "_subblock_vis_stair_case.obj"; // 输出obj的路径

	// 调用函数
	vector<Point3ff> outputVertices;
	vector<face> outputFace;

	// 新版读入
	ifstream ifile(input_file_name);
	int num_contours;
	ifile >> num_contours;
	vector<Point3ff> contour1, contour2;
	for (int i = 0; i < num_contours; i++)
	{
		// 每个contour,读入是否为环isLoop
		int isLoop = -1;
		ifile >> isLoop;
		contour1.clear();
		contour2.clear();
		// 外层contour
		int num_points;
		ifile >> num_points;
		for (int j = 0; j < num_points; j++)
		{
			// cout << "a";
			Point3ff currentP;
			ifile >> currentP.x >> currentP.y >> currentP.z;
			contour1.push_back(currentP);
		}
		if (isLoop)
		{
			// 内层contour
			ifile >> num_points;
			for (int j = 0; j < num_points; j++)
			{
				// cout << "a";
				Point3ff currentP;
				ifile >> currentP.x >> currentP.y >> currentP.z;
				contour2.push_back(currentP);
			}
		}
		printf("正在处理第%d/%d个轮廓或环.\n", i + 1, num_contours);
		newGenerateMesh(vn, isLoop, contour1, contour2, h, outputVertices, outputFace);
	}

	// 输出为obj
	outputToObj(outputVertices, outputFace, outputColor, objPath);
	return;
}
